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ABSTRACT 

An  overview  and  investigation  of  the  more  popular  digital  filter  design  tech¬ 
niques  are  presented,  with  the  intent  of  providing  the  filter  design  engineer  a  com¬ 
plete  and  concise  source  of  information.  Advantages  and  disadvantages  of  the 
various  techniques  are  discussed,  and  extensive  design  examples  used  to  illustrate 
their  application  to  specific  design  problems.  Both  HR  (Butterworth,  Chebyshev 
and  elliptic),  and  FIR  (Fourier  coefficient  design,  windows  and  frequency  sampling) 
design  methods  are  featured,  as  well  as  the  Optimum  FIR  Filter  Design  Program 
of  Parks  and  McClellan,  and  the  Minimum  p  -  Error  HR  Filter  Design  Method  of 
Deczky. 
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Digital  filter  design  today  is  a  rapidly  growing  field  with  myriad  applications 


in  the  areas  of  signal  processing,  image  processing,  filtering,  prediction,  and  estima¬ 


tion.  Extensive  research  has  resulted  in  a  plethora  of  available  information,  which 


can  be  overwhelming  to  the  engineer  presented  with  a  specific  design  problem.  As 


the  title  suggests,  the  purpose  of  this  thesis  is  to  present  an  overview  of  the  more 


popular  design  techniques  in  an  attempt  to  consolidate  the  information  available  in 


the  literature,  and  provide  an  easily  readable  reference  manual.  Detailed  examples 


are  given  to  illustrate  the  application  of  selected  techniques  to  specific  filter  design 


problems. 


Due  to  the  wealth  of  information  available,  a  summary  and  investigation  of 


every  filter  design  technique  is  not  possible.  As  stated  the  more  popular  meth¬ 


ods  are  outlined  in  detail,  however,  references  are  included  as  sources  of  further 


information. 


The  survey  is  threefold,  consisting  first  of  recursive  HR  design  methods,  fol¬ 
lowed  by  nonrecursive  FIR  design,  and  finally  computer-aided  design  (CAD).  A 


summary  of  the  chapter  contents  follows  to  enable  the  reader  to  immediately  lo¬ 


cate  the  section  applicable  to  his/her  particular  design  problem. 


Chapter  II,  Recursive  Filter  Design,  presents  Butterworth,  Chebyshev  and 
elliptic  filter  design  from  both  a  traditional  approach,  wherein  analog  prototype 
filters  based  on  design  specifications  are  converted  to  digital  versions  using  the 
bilinear  transformation,  and  a  direct  design  approach  that  eliminates  the  need  to 


determine  an  analog  prototype. 


Lt'i 


Chapter  III,  Nonrecursive  Filter  Design,  includes  an  overview  of  Fourier  co¬ 
efficient  design,  windowing,  and  a  derivation  and  illustration  of  the  method  of 
frequency  sampling.  Among  the  design  examples  presented  is  the  design  of  a  band¬ 
pass  differentiator  and  a  bandpass  integrator. 

In  Chapter  IV,  Computer-Aided  Design,  the  Optimum  FIR  Filter  Design 
Method  of  Parks  and  McClellan  (that  employs  the  Remez  Exchange  Algorithm) 
is  investigated  and  applied  to  a  bandpass  filter  design  problem.  Also  included  is 
a  short  discussion  of  a  method  to  enhance  the  application  of  this  program  to  the 
design  of  high-order  filters. 

To  complete  the  chapter,  the  Minimum  p- Error  Design  Method  for  HR  Filter 
Design,  that  utilizes  the  Fletcher- Powell  algorithm,  is  presented.  Furthermore,  a 
discussion  of  the  advantages  and  disadvantages  of  these  two  iterative  techniques  is 
included. 


A.  INTRODUCTION 


The  recursive  realization  of  digital  filters  is  advantageous,  in  that  the  desired 
frequency  response  can  be  obtained  using  a  lower  order  filter  than  if  a  nonrecursive 
realization  were  used  (assuming  linear  phase  is  not  required).  This  is  because  the 
filter  frequency  response  is  influenced  by  both  the  poles  and  zeros  of  the  filter 
transfer  function,  whereas  the  frequency  response  of  a  nonrecursive  realization  is 
determined  only  by  the  filter’s  zeros. 

Traditionally,  the  design  of  recursive  digital  filters  involves  the  determination 
of  an  analog  prototype  using  one  of  several  methods:  Butterworth,  Chebyshev  or 
elliptic,  to  name  a  few,  and  converting  these  to  a  digital  version  using  impulse- 
invariant  design  or  bilinear  transformation  methods,  the  latter  being  the  more 
popular. 

This  chapter  first  presents  a  review  of  the  traditional  design  techniques  us¬ 
ing  examples  involving  Butterworth,  Chebyshev  and  elliptic  filters;  then  a  Direct 
Design  method  is  presented,  whereby  the  requirement  for  determining  an  analog 
prototype  is  eliminated,  thus  reducing  the  algebraic  complexity  of  recursive  filter 
design. 

B.  TRADITIONAL  DESIGN  TECHNIQUES 

Recursive  filter  design  techniques,  as  stated  in  the  introduction,  involve  deter¬ 
mining  an  analog  prototype  filter  from  given  filter  specifications,  and  then  convert¬ 
ing  the  prototype  to  a  digital  version  using  a  bilinear  transformation.  Since  this  is 
a  well  known  design  procedure  [1],  details  of  the  derivation  will  not  be  considered. 


c 


However,  as  a  quick  review,  three  examples  involving  the  design  of  Butterworth, 
Chebyshev  and  elliptic  filters  will  be  given. 

It  should  be  noted  that  the  tables  of  analog  prototype  transfer  functions  found 
in  this  section  (Tables  2.2,  2.3  and  2.4),  all  have  a  critical  frequency  of  unity 
( wc  =  1).  Thus,  to  obtain  a  filter  transfer  function  with  a  different  critical  fre¬ 
quency  based  on  these  prototypes  requires  the  use  of  an  appropriate  analog  fre¬ 
quency  transformation  (Table  2.1). 

For  example,  suppose  a  lowpass  analog  filter  with  a  critical  frequency  of  wc 
is  desired.  The  following  relationship  applies: 


HlpUwc)  -  HLPp(jl) 


(2.1) 


where  Hlpp  is  the  prototype  filter  with  a  critical  frequency  of  wc  =  1. 

To  convert  the  prototype  transfer  function  to  the  desired  transfer  function, 
Hlp(s),  the  s  in  the  lowpass  prototype  filter  is  replaced  by  s/wc  as: 


(2.2) 

Table  2.1  (Reference  1)  lists  the  appropriate  frequency  transformations  for  high- 
pass,  bandpass  and  bandstop  filters.  A  summary  of  the  design  procedure  follows. 
1.  Bilinear  Transformation  Design 

From  given  analog  filter  specifications,  determine  the  approppriate  s- 
domain  design  technique  (i.e.,  Butterworth,  Chebyshev  or  elliptic). 

Design  the  analog  prototype  filter  in  accordance  with  the  technique  se¬ 
lected  above.  This  involves  the  translation  of  the  filter  specifications  to  those  of 
a  lowpass  prototype.  The  lowpass  prototype  is  then  designed  according  to  the 
translated  specifications. 


HLP(s)  =  HLPp(s) 
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Transform  the  lowpass  prototype  to  the  desired  lowpass,  highpass,  band¬ 
pass,  or  bandstop  filter  according  to  Table  2.1. 

Transform  the  analog  filter  transfer  function  to  a  digital  version  using  the 


bilinear  transformation. 


H{z)  =  H(s) 


Example  2.1  Butterworth  Filter 

Design  a  digital  filter  that  is  flat  in  the  passband  from  0  to  the  3  dB  cutoff 
frequency,  fc ,  of  2  kHz.  For  frequencies  greater  than  4  kHz,  the  attenuation  should 
be  at  least  10  dB  The  sampling  rate  is  20  kHz. 

Step  1:  A  Butterworth  design  is  called  for  because  a  flat  passband  is  desired. 
Step  2:  Convert  the  critical  analog  design  frequencies  to  digital. 


Sampling  frequency:  fs  =  20  kHz 
Cutoff  frequency:  fc  =  2  kHz 

Stopband  frequency:  /a  =  4  kHz 


w3  =  40  x  1037r  rad /s 
wc  =  4  x  1037t  rad /s 
wa  =  8  x  1037r  rad/s 


ffc  =wcT  =  wjf ,  =  =  °’27r  rad 

&a  =  waT  =  wa/f3  =  ^O^IO^  =  0,47r  rad 
Step  3:  Prewarp  the  analog  frequencies  to  yield  the  desired  digital  frequencies. 

w'c  =  tan  {&c/ 2)  =  tan(0.l7r)  =  0.325  rad/s 

w'a  =  tan  (0^/2)  =  tan(0.27r)  =  0.726  rad /s 

Note:  Prewarping  is  done  to  ensure  the  analog  filter  will  ultimately 

yield  a  digital  filter  with  the  correct  critical  frequencies.  This  is  best 
visualized  in  the  following  graph  (Figure  2.1)  of  the  relationship  between 
the  analog  filter  frequencies  and  the  desired  digital  frequencies, 
i.e.,  w'n  =  tan(0n/2). 


Normalizing  to  the  lowpass  prototype,  according  to  Table  2.1: 

w'  =  0.325  — *  w  =1  rad 
c  Cp 

w'a  =  0.726  — ♦  wap  =  0.726/0.325  =  2.234  rad 


Figure  2.1.  Relationship  Between  Analog  and 
Digital  Filter  Frequencies 


Note:  Normalizing  is  done  to  transform  the  prewarped  critical  frequencies 
of  the  analog  filter  to  their  relative  equivalents  in  a  lowpass  prototype  filter  with 
a  critical  frequency  of  wCp  =  1.  This  is  done  to  enable  the  designer  to  take 
advantage  of  previously  compiled  tables  or  frequency  response  curves  corresponding 
to  the  transfer  functions  of  these  lowpass  prototypes.  Thus,  a  prototype  filter  can 
be  selected  whose  frequency  response  characteristics  have  the  same  shape  as  the 
filter  that  is  being  designed.  The  transfer  function  for  the  selected  prototype  is 
then  converted  to  a  transfer  function  that  exhibits  the  desired  frequency  response 
characteristics  using  an  appropriate  substitution  (Step  5).  Figure  2.2,  depicts  the 
process  of  normalizing. 
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TABLE  2.2 

BUTTERWORTH  PROTOTYPE  COEFFICIENTS 

(Table  after  to  Ref.  4) 


N  01 


1  1.0000 


a2 

03 

04 

05 

06 

1.0000 

2.0000 

1.0000 

3.4142 

2.6131 

1.0000 

5.2361 

5.2361 

3.2361 

1.0000 

7.4641 

9.1416 

7.4641 

3.8637 

1.0000 

10.0978 

14.5918 

14.5918 

10.0978 

4.4940 

13.1371 

21.8462 

25.6884 

21.8462 

13.1371 

Hlpp(s)  = 


1  +  ai5  +  a2.s2  +  . . .  +  a^V 


iy»yi~ 
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Step  7:  Obtain  a  plot  of  the  digital  filter  frequency  response  to  see  if  the  design 
specifications  have  been  met. 

Figure  2.3  shows  that  the  design  does  indeed  comply  with  the  given  filter 
specifications,  in  that  the  cutoff  frequency  is  0.628  =  0.27T  rad  and  the  stopband 
frequency  is  1.257=0.47r  rad  with  gains  of  —3  dB  and  —14  dB,  respectively. 
Example  2.2  Chebyshev  Filter 

Design  a  digital  bandpass  filter  with  the  following  specifications: 

•  1  dB  ripple  in  the  frequency  band  600  to  900  Hz, 

•  sampling  frequency  of  3  kHz, 

•  maximum  gain  of  —40  dB  for  0  <  /  <  200  Hz. 

Step  1:  Convert  the  critical  analog  design  frequencies,  to  the  corresponding 
digital  frequencies,  0{. 

Sampling  frequency:  /*  =  3  kHz 

Lower  ripple  passband  frequency:  f(  =  600  Hz  =>  W(  —  12007T  rad  /  s 
Upper  ripple  passband  frequency:  /„  =  900  Hz  ==►  wu  =  1800tt  rad  /  s 
Stopband  frequency:  fa  =  200  Hz  =►  wa  =  4007T  rad  js 

e't  =  wtT  =  wt/f3  =  =  0.4tt  =  1.26  rad 

^  =  wuT  =  wu/fs  =  =  O.Qtt  =  1.88  rad 

4007T 

&a  =  W*T  =  wa/fa  =  =  0.1337T  =  1.54  rad 

Ripple  band  center  frequency  :  9'0  =  y/0'(9'u  =  0.49tt  =  0.418  rad 
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Figure  2.3.  Frequency  Response  for  Butterworth 
Lowpass  Filter 


Step  2:  Prewarp  the  analog  frequencies  so  that  the  desired  digital  frequency 

characteristics  will  be  achieved. 

w'(  =  tan(<^/2)  =  tan(1.26/2)  =  0.729  ra d/s 

w'u  =  tan  (0^/2)  =  tan(1.88/2)  =  1.369  rad/s 

w'a  =  tan (0^/2)  =  tan(0. 418/2)  =  0.212  rad/s 

Wq  =  =  y,(0.729)(1.369)  =  1.00  rad/s 

Ripple  band  is:  B'  =  w'u  —  w't  =  1.369  —  0.729  =  0.64  rad /s 
From  Table  2.1  convert  the  bandpass  filter  frequencies  to  the  correspond¬ 
ing  lowpass  prototype  filter  frequencies. 


™LPP  = 


*7  '7 

__  wgp  -  w0* 


'2 

W  BP 

-(1.00)' 

(0.64  )w'BP 

Wqp 

WLPp 

w't 

0.729 

-1 

< 

1.369 

1 

< 

0.212 

-7.04 

w'0 

1.000 

0 

Step  3:  Determine  the  order  of  the  Chebyshev  prototype  filter  that  meets  the 
design  requirements. 

For  a  lowpass  Chebyshev  filter  the  magnitude  -  squared  characteristic  is 
given  by: 

HLPp{jw)  =  l  +  e2Cj/{w)  (2'8) 

where  Cn(w)  is  an  Nth-  order  Chebyshev  polynomial  and  e  indicates  the  degree 
of  ripple. 

In  this  example  the  —40  dB  gain  translates  to  M  —  10-2  for  the  desired 
filter  frequency  response  in  the  stopband,  i.e.,  we  want  \H(jw)\  <  10-2  or 
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\H(jw)\  <  10-4.  The  design  specification  for  1  dB  ripple  in  the  ripple  band 

corresponds  to  a  value  for  e  of  0.2589. 

Making  an  initial  guess  for  the  filter  order  of  N  =  2  yields: 


H2(jw)\  = 


1  + 0.2589  (2u,2  -1)2L.=7.04 
1 


1  +0.2589(99.12  -  l)2 

which  is  not  less  than  10~4. 

Proceeding  further  with  N  =  3  yields: 


=  4.01  x  10' 


H3(jw)\  = 


1  +  0.2589 (4u;3  -  3u;a)2  L,=7. 04 
1 


(2.10) 


1  +  0.2589(1395.65  -  21. 12)2 


=  2.044  x  10“6 


which  is  less  than  10-4 .  Therefore,  the  design  needs  can  be  met  using  a  third-order 
Chebyshev  lowpass  prototype  filter. 

Step  4:  Obtain  the  transfer  function  of  the  third  order  lowpass  prototype  filter 
with  1  dB  ripple.  Looking  at  Table  2.3,  yields  the  following  prototype 


transfer  function: 


HlpM  = 


0.491 

s3  +  0.9S8s2  +  1.238s  +  0.491 


(2.11) 


Since  this  is  an  odd  order  filter,  the  constant  0.491  in  the  numerator  was 
selected  to  make  |jfiT(j0)|  =  1.  Recall  that  for  Chebyshev  filter  design,  the  constant 
K  in  the  numerator  of  the  lowpass  prototype  filter  is  selected  on  the  following  basis: 


for  N  odd;  ff(;0)  =  1 


(2.12) 


for  N  even;  if(y’O)  = 


[1-M2]1 
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TABLE  2.3 

CHEBYSHEV  -  PROTOTYPE  DENOMINATOR  POLYNOMIALS 

(Table  due  to  Reference  4) 

0.5-dB  Ripple  Chebvyshev  Filter  (e  =  0.3493) 


1  /  +  2.863 

2  s1  +  1.425/  +  1.516 

3  i»  +  1.253i»  +  1 .335/  -p  0.716  -  (/  ^  0.626X/2  +  0.626/  +  1.142) 

4  s*  -r  1.197/J  *  l.717/»  +  1.023,  +  0.379  -  {/»  -p  0.351/  +  I.064X/2  +  0.845/  +  0.356) 

5  ,/  -  1.I725I/4  4-  1.9374/1  4-  1.3096/1  +  0.7525/  4-  0.1789 

(/  -  0.3623)!(z  -  0.1120)1  4-  1.0116:1{(/  -r  0.2931)1  4.  0.62521] 

6  /‘  -  1.1592/'  *  2.1718/*  -  1.5898/>  +  1.1719/1  +  0.4324/  +  0.0948 

[(/  -  0.0777)1  _  1 .0085 2]t(/  -  0.2121)1  -p  0.7382:K(/  -  0.2898)2  -  0.2702-1 

7  /i  -r  1.151 2/*  —  2.4126/s  4-  1 .869.4/*  4-  1.6479/3  —  0.7556/2  4-  0.2821/  ■+•  0.0447 

(/  -  0.2562)[(,  -r  0.0570)1  -r  1.0064i][(,  -r  0.1597)2  +  0.8001 1][(/  -r  0.2308)2  f 
0.44791] 

8  /»  -  ). (46I/7  4-  2.6567/'  4.  2.1492/5  -  2.1840/4  4-  J.1486/3  4-  0.5736/2  -r 

0.1525/  *  0.0237 

[(/  -  0.0436)1  *  1 ,0050:j((z  -  0.1242)'-  -p  0.8520=][(/  -t-  0.1859)2  -r  0.56932] 

((/  -r  0.2193)1  *  0.19991] 

9  /*  -  1.1426/'  -  2.9027/7  -r  2.4293/'  -  2.7815/5  +  1.61I4/4  -f  0.9836/3  + 

0.3408/1  *  0.0941,  *  0.0112 

(/  -  0.1984)((/  -p  0.0345)1  4-  1.00402]l(/  -4  0.0992)2  +  0.88292][(/  +  0.1520)2  -r 
0.6553 1](( /  -r  0.1864)1  -p  0.34871] 

10  j  1 0  *  l . 1 40 1  / 9  t  3.1499/'  4-  2.7097 /’  4  3.4409/'  -  2.1442/3  4-  I.5274/4  -r 
0.6270/ 3  -  0.2373/1  _  0.0493/  -  0.0059 

(1/  -  0  0279)1  -  1 .003 3 !][(/  -  0.0810)2  -  0.59051  2U(/  -  0.1261)2  -r  0.71832] 

[{/  -  0.1589)1  -  0.46I2:][(/  -p  0.1761)1  0.1589:] 
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TABLE  2.3  (Continued) 

CHEBYSHEV  -  PROTOTYPE  DENOMINATOR  POLYNOMIALS 

n  1.0-dB  Ripple  Chebyshev  Filter  (e  =  0.5089) 

1  z  1.308 

2  /’  -  0.804/  -  0.637 

3  /’  -  0.738/4  -  1.022/  -  0.327  =  (/  -  0.402X/’  -  0.369/  -  0.886) 

4  /*  -  0.7 1 6/ J  -  1.256/’  -  0.317/  -  0.206  =  </’  -  0.210/  -  0.928K/1  -  0.306/  -  0.221) 

j  si  _  0.706S/4  -  1.499S/1  -  0.6935/4  -  0.4593/  -  0.0817 

Is  -  0.2183)11/  -  0.0675)4  -  0.97354)((/  -  0.1766)4  _  0.60164) 
g  z*  _  0  7012/’  -  I.7459/4  •  0.S670/J  —  0.7715/4  -  0.2103/  -  0.0514 

[(/  -  0.0470):  ..  o.98174]I(/  -  0.1283)4  _  0.71S74;[lz  -  0.1753)4  -  0  26304) 

7  s'  -  0.6979/4  -  1.9935/’  -  1.0392/4  -  1.1444/’  -  0.3S25/4  -  0.16661/  -  0.0204 

( s  -  0.!553)[(z  -  0.0346)4  -  0.98674][{/  -  0.0968):  -  0.79I24])(/  -  0.1399)4  - 
0.4391  =] 

8  /•  -  0.6961/7  -  2.2423/‘  -  1.2117/’  -  1.5796/4  -  0.59S2/1  -  0.3387/4  _ 

0.0729/-  0.0129 

[(/  -  0.0265)4  _  0.989841(1/  -  0.0754)4  -  0.8391  :)[</  -  0.1129)4  -  0.56074) 

[1/  -  0.1332)4  -  0.19694] 

9  /»  -  0.6947/*  -  2.4913/4  -  1.3837/*  -  2.0767/’  -  0.8569/4  -  0.6445/’  + 

0.16644/4  -  0.0544/ -  0.0051 

(/  -  0.12061(1/  -  0.0209)4  -  0.9919’][(z  -  0.0003):  -  0.8723 :)[f/  -  0.0924)4 
-  0.647441(1/  ~  0.113414  -  0.34454] 

10  /io  —  0.6937/ 9  -  2.7406/*  —  1.5557/7  -  2.6363/*  —  1.1585/’  —  1.0389/4  — 

0.3178/’  -  0.1440/4  -  0.0233/  r  0.0032 

[1/  -  0.0170)4  ._  0.9935:!(</  -  0.0767)4  _  0.71 1 34)[(z  -  0.0493):  -  0.8962:1 
(1/  -  0.0967)4  _  0.4567-]((z  -  0.1072)’  —  0.15744) 


2-dB  Ripple  Chebyshev  Filter  (e  =  0.7648) 

* 

1  /  +  1.965 

2  /i  -  1.098/  -  1.103 

3  /’  -  0.988/4  -  1.238/  -  0.491  -  (s  -*■  0.494K/*  -  0.490/  -  0.994) 

4  z4  -  0.953/’  -  1.454/4  0.743/  -  0.276  -  </’  -  0.279/  -  0.987)(/»  -  0.674/  -  0.279) 

j  ,)  /.  0.9368/4  -  1.6888/’  -  0.9744/i  -  0.5805/  0.1228 

(s  -  0.2895)((z  -  0.0895)4  0.9901 :]((/  *  0.2342)’  -  0.61194] 

6  /*  j-  0.9282 /’  -  1.9308/4  -  1.2021/’  -  0.9393/’  -  0.3071/  -  0.0689 

(1/  -  0.0622)4  -  0.9934’][(/  -  0.1699)’  -  0.72724]((/  -  0.2321)’  -  0.2662’] 

7  z7  -  0.9231/*  -  2.1761/’  -  I.42S8/4  -  1.3575/’  -  0.5486/’  -  0.2137/  -  0.0307 

(/  -  0.2054)«z  -  0.0457)4  -  0.9953:]((/  -  0.1281)’  -  0.79824][(/  -  0.1851)’  - 
0.4429’] 

8  /•  -  0.9198/7  -  2.4230/*  -  1.6552/’  -  1.8369/4  -  0.8468/’  -  0.4478/’  -r 

0.1073/  -  0.0172 

((/  -  0.0350)4  _  0.9965’I((/  -  0.0997)=  -  0.8448=][(z  -  0.1492)’  -  0.5644’] 

[(/  -  0.1759)4  -  0.19824] 

9  /»  -  0.9175/*  -  2.6709 j7  -  1.8815/*  -  2.3781/’  -  I.20I6/4  -  0.7863/’ 

0.2442/’  -  0.0706/  -  0.0077 

(/  -  0.1593)((z  -  0.0277)4  _  o.99722][(/  *  0.0797)=  -  0.8769’][(/  -  0.1221)=  - 
0.65094)[fz  -  0.1497)4  -  0.3463=] 

10  /">  -  0.91 59/*  -  2.9195/*  -  2.1079/7  -  2.9815/*  -  1.6830/’  -  J.2445/4  - 

0.4554/J  -  0.1825/4  -  0.0345/  -  0.0043 

[(/  -  0.0224)4  _  0.9978=1(6/  -  0.1013)’  -  0.7143’]((z  -  0.0651)=  -  0.9001=] 

((/  -  0.1277)4  -  0.4586’|[(z  -  0.14I5)»  -  0.1580=] 
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Step  5:  Convert  the  lowpass  filter  prototype  to  a  bandpass  filter  prototype. 
From  Table  2.1: 


Hbp(s)  =  Hlpp(s) 


•*+w2 


0.491 


(oil)  +  0-988  (SSi)  +  1-238  (£&)+ 0.491 

0.129s3 


(2.13) 


s6  +  0.632s5  4-  3.507s4  +  1.394s3  +  3.507s2  +  0.6325  +  1 
Step  6:  Determine  the  digital  version  of  the  transfer  function  using  the  bilinear 
transformation. 


Hbp(z )  =  Hbp(s) 


'  r  +  l 


0.129  (gO 


=n3 


'-i)4 

+i ) 

+  1.238 


(f*f)  +  0.632  (f=f)  +  3.507  (f=f)  +  1.394  (*=»  +  3.507  (f^)‘ 

z  —  T 


z  +  1 


+  0.491 


0.011  (z6  -3z4  +3z2  -  1) 
z6  +  2.153Z4  +  1.786z2  +  0.545 


(2.14) 


Step  7:  Confirm  that  the  filter  design  meets  specifications  by  obtaining  a  fre¬ 
quency  response  plot.  Figure  2.4  indicates  the  specifications  for  a  pass- 
band  of  1.26  to  1.88  rad  with  1  dB  ripple  have  been  met. 


Example  2.3  Elliptic  Filter 

A  lowpass  elliptic  filter  with  the  following  specifications  is  desired: 

•  passband  ripple  of  0.5  dB, 

•  passband  ripple-edge  frequency  of  2  kHz, 

•  stopband  gain  should  be  at  most  —20  dB  for  frequencies  greater  than  4  kHz, 
and 

•  sampling  frequency  is  20  kHz. 
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Figure  2.4.  Frequency  Response  for  Chebyshev 
Bandpass  Filter 
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Step  1:  Convert  the  critical  analog  design  frequencies  to  digital  frequencies. 


Sampling  frequency: 


fa  =  20  kHz 


Passband  ripple-edge  frequency:  /i  =  2  kHz 
Stopband  ripple-edge  frequency:  f2  —  4  kHz 


wi  =  4  x  1037t  rad  /s 
u>2  =  8  x  1037r  rad  /  s 


.  (4  X  IQ3)  7T 

0J  =  wxT  =  tui//,  =  2Q  -  1Q3--  =  0.2tt  =  0.628  rad 

.  (8  x  103)  ir 

$2  =  W2T  =  w2/fs  =  2Q  x  1Q3"  =  0.6tt  =  1.88  rad 

Step  2:  Prewarp  the  analog  frequencies  in  order  to  determine  the  appropriate 
lowpass  prototype  filter. 

w\  —  tan(0j/2)  =  0.325  rad /s 
w'2  =  tan  (02/2)  =  1.376  rad/s 

Step  3:  Determine  the  order  of  the  elliptic  prototype  filter  that  meets  the  design 
specifications.  Find  the  ratio,  R,  of  the  stopband  frequency,  w2,  and  the 
passband  frequency,  w[ . 


R  =  =  sis  =  4234 


(2.15) 


From  Table  2.4b  determine  the  prototype  filter  order  and  the  resulting 
transfer  function.  Looking  at  Table  2.4,  it  can  be  seen  that  a  filter  of 
order  N  =  2  has  a  value  for  R  of  2.76261,  which  is  more  than  sufficient 
to  meet  the  design  specifications. 

Step  4:  Obtain  the  transfer  function  of  the  second-order  lowpass  prototype  filter. 

Again,  looking  at  Table  2.4,  the  prototype  transfer  function  is  determined 


to  be: 


HLpP{s)  = 


0.1s2  +  0.5338 
s2  +  0.8094s  +  0.5667 


(2.16) 
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TABLE  2.4a 


ELLIPTIC  FILTERS 
Generalized  Transfer  Functions 

(due  to  reference  5) 


I 

i 

> 

i 

r 
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i 
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} 
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N  =  2 


N  =— 3 


N  =  4 


TT  (  \  _  T r  -S2  +  Aqi _ B032  +  Up  Aqi 

32  +  B\\a  4-  Bq\  s 2  +  Bn  3  4-  Roi 

TT  /  X  _ _  Bq  S  Ap  1 

3  +  So  32  +  Bn  3  +  Boi 

JJ  /  X _ HqS2  4~  HqAqi _ 

a3  4-  (Bn  +  3p)  ^2  +  (Boi  4-  B\\ sp)  s  +  Bpi  sp 
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JJ  /X  _  _ HqS4  -f  Hq  (Aqi  +  A02)  •S2BpApx  A02 _ 
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4*  [BxiBp2  4-  BpxBiz  4-  B02S0  4-  B\iB\2Sq  4-  BoxSo]  s2 

4-  [B01 B02  4-  B\\ Bp23p  4-  Bpx  Bi23p]  s  4-  Bpx  Bp23p 
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ELLIPTIC  FILTERS 
(due  to  reference  5) 
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TABLE  2.4c 

ELLIPTIC  FILTERS 

(due  to  reference  5) 
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TABLE  2.4d 

ELLIPTIC  FILTERS 

(due  to  reference  5) 
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Step  5:  Find  the  analog  transfer  function  of  the  desired  filter  based  on  the  lowpass 
prototype. 

This  is  done  by  finding  a  scaling  factor,  a,  to  convert  the  lowpass  proto¬ 
type’s  cutoff  frequency,  uqp,  to  the  desired  cutoff  frequency,  w\ . 

For  the  prototype: 

R  =  2.76261  ;  wlp  =  1/VR  =  0.6016 


For  the  desired  filter: 


to  i  =  0.325 


Therefore,  the  scaling  factor,  a,  is: 


/  0-6016  1 
a  =  u,1'/u,1  =  a^  =  L85 


(2.17) 


Thus,  the  actual  transfer  function,  Hip(s),  can  be  found  as  follows: 
BLF(s)=HlPr(s)  =HLPf(s) 

s—cx  s  j=  1 .8  5« 

0.1s2  +0.5338 

s2  +  0.8094s  +  0.5667  a=185a  (2.19) 

0.1(1.85s)2  +0.5338 
“  (1.85s)2  +0.S094(1.85s)  + 0.5667 
0.3423s2  +  0.5338 
“  3.423s2  +  1.497s +  0.5667 

Step  6:  Determine  the  digital  version  of  the  transfer  function  using  bilinear  trans¬ 
formation. 


Hlp(z)  =  Hlp(s) 


0.3423  (z  -  l/z  +  1)  +  0.5338  _ 


3.423 (z  -  1/2  +  l)2  +  1.497(2  -  1  /z  +  1)  +  0.5667  (2.20) 
0.8761z2  +  0.3832  +  0.8761 
5.487z2  —  5.7126s  +  2.493 
0.1597z2  +  0.0698z  +  0.1597 
z2  -  1.0411z +  0.4543 
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Step  7:  Verify  that  the  design  meets  the  specifications  from  a  frequency  response 
plot  (Figure  2.5). 

2.  Summary 

This  concludes  the  review  of  the  use  of  the  bilinear  transformation  to 
design  recursive  filters  based  on  Butterworth,  Chebyshev  or  elliptic  analog  pro¬ 
totypes.  As  can  be  seen,  this  method  is  very  involved  in  terms  of  algebraic  ma¬ 
nipulation.  The  following  sections  will  introduce  the  direct  design  technique.  It 
will  be  shown  that  the  direct  design  procedure  reduces  considerably  the  amount  of 
algebraic  calculations  required,  and  eliminates  the  somewhat  confusing  procedure 
of  prewarping. 


C.  DIRECT  DESIGN  TECHNIQUE 


The  direct  design  technique  is  based  on  the  bilinear  transformation  in  the 
following  manner: 

The  tables  of  coefficients  for  analog  lowpass  prototype  filters,  i.e.,  Butterworth  - 
Table  2.2,  Chebyshev  -  Table  2.3,  and  elliptic  -Table  2.4,  were  converted  to  z- 
domain  versions  through  the  use  of  the  bilinear  transformation,  and  can  be  shown 
to  have  a  critical  frequency  of  7r/2.  The  reason,  as  stated  in  Section  A,  is  that 
the  analog  prototype  filters  comprising  these  tables  all  have  a  critical  frequency  of 
wc  =  1.  For  the  bilinear  transformation,  the  following  substitution  for  s  is  used  in 
the  analog  prototype  transfer  function  to  convert  it  to  a  digital  prototype  transfer 
function: 


s 


z  —  1 
2  +  1 


or  z  = 


1  +  s 
1  —  s 


(2.21) 


An  analog  critical  frequency  of  s  =  jwc  =  j  1  implies: 


o.ooooo  o.r.2nn2 


2.51326 


3.14160 


1.25G64  1.88-196 

FREQUENCY,  THETA 


N  =  2 


Figure  2.5.  Frequency  Response  for  Elliptic 
Lowpass  Filter 
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since  z  =  ,  this  gives: 


» Jec  =  — —  =  \ei*l2 

1  ~  j 


(2.23) 


Thus,  it  can  be  seen  that  an  analog  critical  frequency  of  wc  =  1  yields  a 
digital  critical  frequency  of  dc  —  tt/2.  Tables  of  prototype  filters  were  thus  created, 
specifically,  Butterworth  -  Table  2.7,  Chebyshev  -  Table  2.8  and  elliptic  -  Table  2.9. 

The  use  of  these  tables  to  design  a  digital  filter  based  on  analog  specifications 
involves  the  following  generalized  steps.  To  further  illustrate  the  technique,  a 
detailed  description  of  its  derivation  and  use,  including  examples  for  the  various 
filte)  types,  will  be  presented. 

1.  Direct  Design 

Step  1:  Determine  the  appropriate  desired  filter  type  based  on  the  design  speci¬ 
fications,  i.e.,  Butterworth,  Chebyshev  or  elliptic. 

Step  2:  Convert  the  analog  filter  frequency  specifications  to  digital  equivalents. 
Step  3:  Use  the  digital-digital  frequency  transformations  of  Table  2.5  to  normal¬ 
ize  the  desired  filter’s  design  frequencies,  so  that  the  appropriate  lowpass 
prototype  filter  may  be  selected. 

Step  4:  From  the  normalized  design  frequencies  obtained  in  Step  3,  determine 
the  lowpass  prototype  that  meets  or  exceeds  these  requirements. 

For  Butterworth  and  Chebyshev  filters  -  this  is  done  using  the  design 
curves  of  Figures  2.8  and  Figures  2.10  -  2.12. 

For  elliptic  filters  -  Table  2.9  is  used. 

Step  _5:  Obtain  the  actual  filter  transfer  function  from  the  lowpass  prototype, 
based  on  the  digital-digital  frequency  transformations  of  Table  2.5. 
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TABLE  2.5 

DIGITAL  FILTER  FREQUENCY  TRANSFORMATIONS 

(Due  to  Reference  2) 


Type  Transformation 

(Replace  z  in  LP  digital 
prototype  with) 


1.  Lowpass 


1  — a* 


Design  Constants 


_  «in(gc/2-^/2) 
Sin(0c/2+*'c/2) 


2.  Highpass 


1  —  az 


cos(gc/2-^/2) 

cos(*c/2+0'c/2) 


3.  Bandpass 


_  co9(0>M/2+tf'</2) 

“  cos(8>j2-9't/2) 

k  =  tan  ^p-cot  {0'u/2  —  0'(/ 2) 


4.  Bandstop 


«3-.T»fr«+fefc  =  co*(K/2+W) 

cos(^/2-<^/2) 

fc  =  tan  tan  (0^/2  —  0^/2) 


V  v. 


As  can  be  seen  from  the  previous  design  summary,  the  crux  of  this  method 
is  the  use  of  the  digital-digital  frequency  transformations  of  Table  2.5,  which  merit 
explanation  [2]. 

These  transformations  enable  the  user  to  convert  the  lowpass  digital  pro¬ 
totype,  HlPp(z)  (be  it  Butterworth,  Chebyshev  or  elliptic)  to  the  actual  lowpass 
( LP ),  highpass  (HP),  bandpass  (BP),  or  bandstop  (BS)  filter  transfer  function. 

The  transformations  provide  a  means  of  transferring  the  stability,  inherent 
in  the  prototype  filter,  to  the  actual  filter;  that  is,  the  poles  of  the  actual  filter  will 
lie  inside  the  unit  circle,  as  they  do  for  the  prototype  filter.  For  this  reason  the 
frequency  response  of  the  lowpass  prototype  filter  at  a  specific  frequency  of  9a  must 
be  the  same  as  the  desired  value  of  the  frequency  response  of  the  actual  filter  at 


its  corresponding  frequency  of  9la. 


Hlpp  =  Ht„,  (e*) 


(2.24) 


where, 


9a  =  prototype  frequency 
8'a  =  desired  filter  frequency 
Type  =  LP,  HP,  BP,  or  BS 

For  example,  when  transforming  a  lowpass  prototype  filter  to  a  desired 
lowpass  filter,  one  wishes  to  maintain  the  integrity  of  the  magnitude  characteristic 
of  the  prototype  filter,  while  expanding  or  compressing  the  frequency  scale  so  that 
the  critical  frequency  is  changed  from  the  prototype  value  of  9C  =  n/2  to  the 
critical  frequency  of  the  actual  filter  9'c,  (see  Figure  2.6). 

According  to  Table  2.5,  this  involves  the  following  substitution  for  z  in 


the  lowpass  prototype  transfer  function. 


1  —  q: 


(2.25) 


Figure  2.6.  Lowpass  Prototype  -  Lowpass  Transformation 
For  the  critical  frequency  of  the  prototype,  9C ,  to  map  to  the  critical 
frequency  of  the  actual  filter  ,  the  following  must  be  true: 

Substituting,  =  z  on  the  left  in  Equation  (2.25),  and  e*$c  —  z  on  the 


right  gives, 


ej9c  = 


—  a 


1  —  ae}Sc 

and  solving  for  the  design  constant,  a,  yields, 


sin  ( 9c/2-9'r/2 ) 
sin(9c/2  +  9'c/2)- 


Similar  arguments  apply  to  HP,  BP,  and  BS  filters. 


(2.26) 


Thus,  the  transformations  of  Table  2.5  accomplish  the  appropriate  map¬ 
ping  from  the  lowpass  prototype  digital  filter  to  the  actual  filter,  and  ensure  that 
stabilty  is  maintained. 
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2.  Butterworth  Filters  -  A  Direct  Design  Approach 

The  design  curves  shown  in  Figures  2.8a  -  d,  describe  lowpass  prototype 
filters  with  a  critical  frequency  of  9C  =  ir/2,  and  of  order  N  =  1  through  N  =  8. 
Higher  order  filters  cam  be  realized  by  cascading  the  appropriate  lower  order  filters 
characterized  by  these  curves. 

The  curves  were  constructed  by  using  the  s-domain  version  of  the  But¬ 
terworth  filter  coefficients  for  N  =  1  through  N  =  8,  and  then  determining  the 
^-domain  version  of  these  transfer  functions  through  the  use  of  the  bilinear  trans¬ 
formation  and  the  conversion  tables  obtained  from  Reference  3,  (see  Table  2.6). 
These  tables  establish  the  necessary  relationships  between  the  coefficients  of  an 
analog  filter  and  its  digital  counterpart. 

a.  Generalized  Analog  Transfer  Function 


H(s)  = 


Ao  +  Ais  +  A2S2  +  . . .  +  Aksk 
Bq  +  B\S  +  B2&2  +  . . .  +  BkSk 


(2.27) 


b.  Generalized  Digital  Transfer  Function 

_  a0  +  aiz-1  +  Q2 z~2  +  .  ■  •  + 

1  +  biz~l  +b2z~2  +  . . .  +  bkz~k 


(2.28) 


For  example,  in  the  case  of  N  =  1  the  bilinear  transformation  digital 
filter  coefficients,  in  terms  of  the  analog  filter  coefficients  are: 

A  Bq  +  B\ 
ao  (-^o  +  A\ )/ A 

a  1  (Ao  —  A\)/A 

b 1  (B0-B1)/A 

Note:  C  is  the  critical  frequency,  wc  =  1. 

For  a  first-order  Butterworth  prototype  filter  described  by: 
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TABLE  2.6 

BILINEAR  TRANSFORMATION  DIGITAL  FILTER  COEFFICIENTS 
IN  TERMS  OF  ANALOG  FILTER  COEFFICIENTS 


A  8g*B1C 
*0  (Ao+*iCl/A 
o,  (a^CI/A 
b,  <80-8,0/* 

(1st  order) 


a  Bo-a,c»92c2*8jC3-84c4 

Qq  (Aq.-A.  C*  A2C2-AjC3*A4C41 /A 
0,  (AA0-2A,C-2AjC3-AA4C4)/A 

o2  <6A0-2A2C2-6A4C4)/A 
Bj  <4A0-2A,C*2AjC3-AA4C4:/A 

a4  <a0-a,c+a2c2-a3c3*a4C4!/a 
t>,  (ABo-28.  C-293C3-AB4C4! /A 

6,  <63o-2B2C2*-684C4)/a 

bj  (ABo-28,  C»2B3C3-J84CJ) /A 
b4  (Bo-a^^Bz^-Sj^^C'l/A 

(4th  order) 


(due  to  reference  3) 

A 

Bo»B, C»8jC2»8jC3 

b0*b1c+82c2 

°° 

(A0*A,  C«-A2C2't-A3C3)  /  A 

<a0*a,c*a2c2i/a 

(JAq-^C-A^-JAjC^/A 

(2A0-2A2C2)/A 

°2 

<3A0-A,C-A2C2*3AjC3)/A 

(A0-A,C»A2Cil/A 

°3 

(Ao-A, C*A2C2-A jC1 )  /A 

<2B0-2B2C2)/A 

b1 

(3B0*81C-B2C2-3B,C3)/A 

(ao-8,C*B2C2)/A 

(380-8, C-B2C2»383C3)/A 

(2nd  order). 

bJ 

(B0-B1C't-B2C2-B3C3)/A 

(3rd  order). 

A  Bo*B,C*B2C**05C3*B4C4*B,CS 
a0  [  A0*A,  C*A2C*«-AjC3  «-A4C4*A,CS )  /A 

а,  (5A0«JA.C+A2C2-AjC3-1A4C4-5A5C3)/A 

a2  <  lOAo*2A,  C-2A2C?-2AjC3+2A4C4  *  iOA,C3)  /A 
as  ( I0A0-2A, C-2A?C2 ♦2a3C3  +2a4C4 - iOA,c5)/A 
o4  <5A0-JA,C‘A2C?*AjC3-3A4c4*5A5C5)/A 
as  <A0-A,CfA2C2-AjC3*A4C4-AjC3)/A 

б,  <5Bo»3B,C+02C2-BjC3-SB4C4-595C3)/A 

b2  <lOBo*2B,C-2e2C2-2BjC3*2B4C4*tOB5CS)/A 
bj  <tO80-2B,C-2B2C2*2BjC3+2B4C4-tOB5Csl/A 
b4  <5B0-iB,C+B2Cz+B3C3-3B4C4*585C3)/A 
bj  <B0-B,C^B2C?-BjC3  +  B4C4-B5C3)/A 

(5  th  order) 


the  values  for  the  digital  filter  coefficients  axe: 

A  1+1=2 

a0  (l+0)/2  =  0.5 

aa  (1  -  0)/2  =  0.5 

bi  (1  -  l)/2  =  0 
The  transfer  function  H(z)  is  therefore: 

#i(z)  = - - - =►  Ex{z)  = -  (2.30) 

1  z 

A  plot  of  this  transfer  function  is  illustrated  by  the  curve  (N  =  1)  in  Figure  2.7 


Table  2.7  is  a  summary  of  the  analog  and  digital  versions  of  lowpass  prototype 
Butterworth  filter  transfer  functions  for  filter  orders  up  to  N  =  8. 

To  illustrate  the  direct  design  approach  the  design  problem  of  the 
previous  section  will  be  used  (Example  2.1),  where  the  problem  statement  was: 

Design  a  digital  filter  for  a  20  kHz  sampling  rate  that  is  flat  in  the 
passband  of  0  to  the  3  dB  cutoff  frequency  of  2  kHz,  and  has  a  gain  of  not  more 
than  —10  dB  for  frequencies  greater  than  4  kHz. 

Step  1:  A  Butterworth  design  is  called  for  because  a  flat  passband  is  desired. 
Step  2:  Convert  the  critical  analog  design  frequencies  to  digital. 

Sampling  frequency:  f3  =  20  kHz 

Cutoff  frequency:  fc=2  kHz  =>  wc  =  4  x  1037r  rad  j  s 

Stopband  frequency:  fa  =  4  kHz  =►  wa  —  8  x  1037r  rad/s 

6'  =u,cT  =  u,c//.  =  liii^=0.2,rrad 

c  c  c /  20  x  103 

.  (8  X  103)  7T 

#a  =  waT  =  wa/f.  =  2Q  x  103  =  0.47T  rad 

Step  3:  Since  the  “design  chart”  is  based  on  a  critical  frequency  of  0C  —  7t/2, 
digitized  design  specifications  need  to  be  normalized  to  this  frequency, 
so  that  the  order  of  the  required  lowpass  prototype  filter  can  be  deter¬ 
mined.  For  this  process  the  following  frequency  transformations  apply 
(from  Table  2.5,  prime  denotes  desired  filter). 

•  Prototype  filter  from  desired  filter: 


TABLE  2.7 

LOWPASS  PROTOTYPE  BUTTERWORTH  FILTER 
TRANSFER  FUNCTIONS 
(due  to  reference  4) 


Filter  Order 


HlpM 


a2  +  1.4142a+l 


a3+2a2+2a+l 


a4 +2 .6131  a3 +  3.4 14  2a2  +2.6131  a+1 


a5 +3 .23  61  a4 +5. 2  36  la3 


+5 .2  36  la2 +3. 23  61  a+1 


a® +3 .8637a® +7 .4641  a4 +  9.1 41 6a3 +  7.4641  a2 


+3. 8637  a+1 


a7 +4 .49  40  a®  +  10.0978as+14.5918a4 
+  14.5918a3  +  10.0978a2+4.4940a+l 


a8 +5. 1258  a7  +  l  3. 131 7a® +21. 8462  as 
+  25.6884a4+21.8462a3  +  13.1371a2 


+5.1258a+l 


Hlpp{z) 


0.5(z+l) 

z 

(*+l)2 

3.41z2+0.59 

(*+Q3 

6z3  +  2z 

10.6383z4+5.17z^ 

_ Le±22! _ 

18.94z»  +  11.996z3+1.0549z 


_ (*+ir _ 

33. 79  7z  6 +26. 284  z4 +3. 86  z2 +0.0592 


_ _ Liilil _ _ 

60.367z7+55.537z*  + 11 .633z3+0 .4624  z 


107.896zs  +  l  14. 438  z® +  3 1.4  96  z4 


+2.162zz+0.008 


,  '  y.f. 


Figure  2.8c.  Butterworth  Filter  Design  Curves 
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•  Desired  filter  from  prototype  filter: 


pid*  = 


eJ  »  +  a 
1  + 


where 


(2.32) 


sin  (Sc/2-«^/2) 

sm(ec/2  +  6'c/2) 


ffc  and  ffa  are  the  critical  and  stopband  frequencies  (respectively)  used  to 
determine  the  optimum  prototype  filter  order. 

Equation  (2.31)  is  applicable  to  this  problem,  and  thus: 

da  =? 

9'a  =  0.47T 
9C  =  0.57T 
=  0.2ir 

_  sin(Q.57r/2  —  0.27t/2) 
a  sin(0.57r/2  +  0.27r/2) 

_  sin(0.157r)  _  0.454  = 
sin(0.357r)  0.891 


Therefore, 


e>°-4,r  -  0.510 
2  *  “  1 -o.5i9e>0-4*- 
0.97ejlJ79 


0.97e">°'57 


=  ej2-349 


(2.33) 


0a  =  2.349  rad  =  0.757T  rad 


Step  4:  Once  the  normalized  stopband  frequency  is  obtained,  (in  this  case  0.757r), 
the  design  chart  is  used  to  determine  the  lowest  order  filter  that  will  give: 


—  3 dB  at  9C  =  0.5?r 


less  than  —  10 dB  at  9a  =  0.7n 
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Turning  to  Figure  2.8,  it  can  be  seen  that  the  curve  for  N  =  2  meets 
these  requirements.  At  this  point,  it  is  worth  noting  that  a  second-order  Butter- 
worth  filter  was  also  determined  to  be  required  using  the  traditional  approach  in 
the  previous  section. 

In  Table  2.7  it  can  be  seen  the  prototype  low-pass  filter  transfer  func¬ 


tion  for  N  =  2  is: 


(z  +  1) 


z2  +  2z  +  1 


HLpP(z)  3.41.j2  _j_  Q  59  3.4l22+0. 


(2.34) 


Step  5:  Once  the  prototype  filter  is  obtained  from  the  design  curves,  it  is  neces¬ 
sary  to  determine  the  actual  transfer  function  that  will  meet  the  design 
specification,  such  that  &  is  0.2 tt  and  ffa  is  0.47T. 

Using  the  digital  filter  frequency  transformations  of  Table  2.5,  the 
transfer  function  is  determined  as  follows: 

Hlp(z)  =  Hlpp(z)  I 


z2  +  2z  +  1 

3.4122  +0.59  2—0.510 

1  —  0 .51  Ox 

~  *2  +  2z  +  1 

~  14.8422  -  17.02  +  6.14 

which  is  the  same  transfer  function  obtained  in  the  previous  section. 


(2.35) 


Step  6:  A  computer  generated  frequency  response  plot  of  the  filter  is  used  to 
verify  that  the  design  fulfills  the  design  specifications;  Figure  2.9  reveals 
that  this  is  indeed  the  case. 

D.  LOWPASS  PROTOTYPE  TO  HIGHPASS,  BANDPASS  BAND- 
STOP  FILTERS 

Once  the  lowpass  prototype  filter  is  obtained,  a  highpass,  bandpass  or  band- 
stop  filter  can  be  derived  using  the  digital  frequency  transformations  of  Table  2.5. 
In  the  next  section,  an  example  is  presented  that  illustrates  this  process. 
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1.  Chebvshev  Filters  -  A  Direct  Design  Approach 


As  stated  in  the  introduction  to  this  section,  this  design  approach  is  of 
course  applicable  to  Chebyshev  filter  design  problems.  Again  using  the  bilinear 
transformation  techniques  cited  in  [2],  the  digital  filter  coefficients  corresponding 
to  the  analog  filter  coefficients  for  prototype  Chebyshev  filters  were  determined. 
The  frequency  responses  of  the  resulting  transfer  functions  listed  in  Table  2.8  were 
then  plotted  to  obtain  a  set  of  design  curves  with  a  cutoff  frequency  of  7r/2  (see 
Figures  2.10  thru  2.12). 

The  design  problem  of  the  previous  section  will  be  used  to  illustrate  ap¬ 
plication  of  this  new  technique  (Example  2.2). 

The  problem  statement  called  for  a  Chebyshev  digital  bandpass  filter  to 
be  designed  to  meet  the  following  specification: 

•  1  dB  ripple  in  the  range  of  600  to  900  Hz 

•  Sampling  frequency  is  3  kHz 

•  Maximum  gain  of  — 40dB  for  0  <  /  <  200  Hz 

Step  1:  Convert  the  critical  analog  design  frequencies  to  digital: 


Sampling  frequency: 


fs  =  3  kHz 


Lower  Passband  frequency:  ft  =  600  Hz 
Upper  Passband  frequency:  /„  =  900  Hz 


Stopband  frequency: 


fa  =  200Hz 


*t  =  vtT  =  wt/fs  =  =  0.47T  =  1.26  rad 

0'u  =  wuT  =  wu/fs  =  3QQQ  ~  =  0,67r  =  li8S  rad 
9'a  —  waT  =  wa/fs  =  =  0.1337T  =  0.418  rad 


Ripple  band  center  frequency:  0'Q  =  \f0't6'u  =  \/2.35  =  1.54  rad 
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Figure  2.10b.  Chebyshev  Filter  Design  Curves 
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Figure  2.10d.  Chebyshev  Filter  Design  Curves 
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Figure  2.11c.  Chebyshev  Filter  Design  Curves 
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Figure  2.12b.  Chebyshev  Filter  Design  Curves 
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Step  2:  Again,  the  design  chart  is  based  on  a  ripple  edge  frequency  of  tt/2,  necessi¬ 
tating  the  normalization  of  the  design  specification  frequencies.  Looking 
at  the  Table  of  Digital  Filter  Frequency  Transformations  (Table  2.5),  the 
following  frequency  transformations  apply  for  a  bandpass  filter: 


Prototype  filter  from  desired  filter: 


= 


j20'a  _  2ai  ,  k=l 

_  __ _ k± 1  e  Z. 


i  -  wt";  + 


(2.36) 


where: 


_  cos(<?u/2  -f  fft/2) 

“  cos(«i/2-tfJ/2) 

fc  =  tan  (t)c  /2)  cot  7 y  -  -j) 

#a  is  the  stopband  frequency  used  to  determine  the  filter  order  N. 

6C  is  the  ripple  edge  frequency  on  the  design  chart,  in  this  case  9C  =  7t/2. 
For  this  problem: 


6a  =? 


=  0.1337T 


6'u  =  0.6tt 


0't  =  0.4tt 


<?c  =  0.57T 


COS  (0.37T  +  0.27t)  cos(0.57t) 
01  COS  (0.37T  —  0.27t)  cos(0.l7r) 
k  =  tan(7r/4)  cot(0.l7r)  =  3.07 


*  '.ap/  v  v  v 


-  > '  v.V, 


(2.37) 


+  0-509 

“  1  +  ~  ~  1  +  (0.509)e^ 

— gjo.836  _  0  509  1.394e-j2'58 

1  +  0.509e;O-836  “  1.394e>0-282 


=>  9a  =  —2.86  rad  =  2.86  rad  (due  to  symmetry) 

Step  4:  Determine  the  lowest  order  Chebyshev  lowpass  filter  prototype  from  the 
design  chart  that  gives: 

—  3  db  at  0C  =  0.57T  rad 
less  than  —  40  db  at  6a  =  2.86  rad 

Turning  to  Figure  2.11,  the  filter  order  is  determined  to  be  N  =  3.  Again, 
this  agrees  with  the  order  obtained  using  the  approach  in  the  previous  section. 
From  Table  2.8  the  lowpass  prototype  transfer  function  for  N  =  3  is: 


0.49lz3  +  l.473z2  +  1.473*  +  0.491 
LW)  -  3  71?23  _  x.277*2  +  2.247z  -  0.759 


(2.38) 


Step  5:  As  with  the  Butterworth  design,  the  Digital  Filter  Frequency  Transfor¬ 
mation  Table  2.5,  is  used  to  convert  the  lowpass  prototype  filter  to  a 
bandpass  filter. 


Hbp(z) 


HLPp(z ) 


,2_  2cLk.,+  L=± 

i  -  k-l  ,1 

1  k+l  +*+l 


0.491z3  +  1.473z2  +  1.473*  +  0.491 
3.717 z3  -  1.277z2  +  2.247z  -  0.759 
0.0l2z6  -  0.036z4  +  0.036z2  +  0.012 
z6  +  2.136 z4  -I-  1.768z2  +  0.540 


(2.39) 


Again,  this  is  nearly  identical  to  the  transfer  function  obtained  using  the 


“traditional  approach”. 


Step  6:  Verification  of  design  using  a  computer  generated  frequency  response 
plot.  From  Figure  2.13  it  can  be  seen  the  design  specifications  of: 

9\  =  0.4*  =  1.257  rad 
=  0.67T  =  1.885  rad 
e'a  =  0.133tt  =  0.418  rad 

are  met. 

2.  Elliptic  Filters  -  A  Direct  Design  Approach 

As  was  shown  in  Example  2.3,  the  design  of  elliptic  filters  is  algebraically 
intensive,  an  extremely  undesirable  feature,  especially  for  high-order  filters.  Elliptic 
filters  characteristically  exhibit  a  sharper  cutoff,  that  is,  a  narrower  transition  width 
for  a  given  filter  order  than  their  Butterworth  or  Chebyshev  counterparts.  For  this 
reason  they  are  very  popular,  making  the  pursuit  of  a  more  efficacious  design 
approach  a  worthwhile  endeavor. 

Just  such  an  approach  was  found  by  taking  the  analog  prototype  filters  of 
Table  2.4  and  normalizing  them  to  have  a  passband  ripple  edge  (critical)  frequency 
of  one.  They  were  then  transformed  to  digital  versions  using  the  bilinear  transfor¬ 
mation,  resulting  in  digital  prototype  filters  with  a  critical  frequency  of  7r/2  (Table 
2.9). 

Thus,  the  amount  of  calculation  involved  in  the  design  of  a  digital  elliptic 
filter  is  cut  in  half  because  the  requirement  for  finding  an  analog  prototype  is 
eliminated.  The  designer  need  only  convert  the  filter  design  specifications  to  digital 
frequencies,  find  the  digital  lowpass  prototype  filter  that  meets  these  requirements 
(from  Table  2.9),  and  then  find  the  actual  filter  transfer  function  through  the  use 
of  the  digital  frequency  transformations  of  Table  2.5. 
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Two  examples  will  be  presented  that  illustrate  this  process;  however,  a 
more  detailed  explanation  of  how  the  transfer  functions  of  Table  2.9  were  arrived 
at  is  in  order. 

For  purposes  of  explanation,  a  second-order  filter  with  0.5  dB  passband 
ripple,  and  a  stopband  gain  of  —20  dB  will  used. 


Turning  to  Table  2.4,  the  lowpass  prototype  transfer  function  is  of  the 


form: 


where, 


HLpp{s)  = 


HqS2  +  HqAqi 
s2  +  Bus  -f  Boi 


H0  =  0.100220 
Aoi  =  5.33789 
Boi  =  0.566660 
Bu  =  0.809390 


therefore, 


HLPp(s)  = 


Q.ls2  +  0.5338 
s2  +  0.8094s  +  0.5667 


(2.40) 


This  prototype  has  a  value  for  R  of  2.76261.  The  parameter  R  is  the  ratio 
of  the  stopband  frequency,  u;2p,  to  the  passband  cutoff  frequency,  uqp(i.e.,  R  = 
ttf2p  /  w\ p).  Thus,  it  describes  the  sharpness  of  the  transition  region.  A  large  value 
of  R  is  indicative  of  a  broad  transition  region,  while,  conversely,  a  small  value 
corresponds  to  a  narrow  transition  region.  As  expected,  the  higher  the  filter  order, 
the  narrower  the  transition  region,  and  the  smaller  the  value  for  R.  Figure  2.14 
(due  to  [5]),  illustrates  the  magnitude  squared  frequency  response  of  the  normalized 
lowpass  elliptic  filters  of  Table  2.4. 

For  this  particular  example,  wip  and  u>2p  are  subject  to  the  constraint 
that  their  ratio  R  be  2.76261.  Furthermore,  it  should  also  be  noted,  the  filters 
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Figure  2.14.  Magnitude  Squared  Frequency  Response 
of  a  Normalized  Lowpass  Elliptic  Filter 

comprising  this  table  have  been  normalized,  so  that  the  geometric  mean  of  wip 
and  w-ip  is  one. 

(wipw-i,,)  =  1  (2-41) 

Combining  these  two  constraints,  the  following  relationships  between  and  w2p 
apply: 


w\p  =  1  /VR  —  0.6016  rad/s 
w2p  =  y/H  =  1.6620  rad /s 


(2.42) 

(2.43) 


The  goal  of  the  direct  design  procedure  is  to  find  a  digital  parallel  to  the 
analog  prototype  filters,  and  tabulate  a  table  of  digital  lowpass  prototypes. 

Step  1  of  this  process  is  to  normalize  the  transfer  functions  of  Table  2.4 
so  that  they  all  have  a  value  for  u;iP  of  one.  This  is  depicted  Figure  2.15. 

For  this  example: 

0.1s2  +0.5338 

Lpp(3)  s2  +  0.8094s  +  0.5667  ^  ^ 


(2.44) 


Original 


Normalized 


where 


Figure  2.15.  Normalization  of  Elliptic  Filter 


wlp  =  l/y/R  =  1/^2.76261  =  0.6061 


To  normalize  Hlpp($ )  to  w\p  =  1,  let  3  =  s/\/R  =  0.60615 

0.1(0.606l5)2  +  0.5338 

LPp  *  ”  (O.6O6I5)2  +  0.8094(0.60615)  +  0.5667 
0.036752  +  0.5338 
~  0.367452  +  0.49065  +  0.5667 


(2.45) 


where  now,  wlp  =  1.  and  tv2p  =  R,  since  R  =  w^p/w \p. 

The  normalized  analog  prototype  was  then  converted  to  a  digital  prototype 
through  the  use  of  bilinear  transformation. 


0.0367  (f^j)  +0.5338 


0.367  (f^[)2  +  0.4906  (^j— )  +  0.5667 


(2.46) 


1.575Z2  -  2.75*  +  1.575 
3.91*2  +  1.13*  +  1.22 

The  digital  passband  critical  frequency  is  9Xp  =  m/2,  since  an  analog 


frequency  of  uqp  =  1  corresponds  to  digital  frequency  of  tt/2  when  using  the 


bilinear  transformation. 
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Table  2.9  contains  the  compiled  results  of  applying  this  process  to  some  of 
the  analog  prototype  transfer  functions  of  Table  2.4.  As  with  analog  elliptic  filters, 
special  relationships  exist  between  the  passband  critical  frequency,  9lp,  and  the 
stopband  frequency,  02p,  of  digital  elliptic  filters. 

Using  the  digital  frequency  transformations  of  Table  2.5,  the  following 
values  for  0Jp  and  62 P  are  found  based  on  the  analog  prototype  values  for  w\p  and 
w2p. 


For  0i  p, 


=  Jl± 1 


le^/4 
-j  1  +  1  ~  le“-?ir/4 


=  lej*/2 


(2.47) 


For  02  P , 


01 P  =  tt/2 


,je2p  _  iff  +  1 


le* 


I  tan  1  R 


-jR+  1  1  e-;tan-*(-/t) 

02  p  =  2  tan-1  R 


=  lej2 


tan”1  R 


(2.48) 


Again,  R  is  the  transition  region  parameter  described  by  the  ratio  of  w2p 
to  W\p  of  the  analog  version  of  the  filter.  Table  2.9  presents  a  summary  of  the 
values  for  0 ip,  02p ,  and  R  for  some  of  the  digital  lowpass  prototype  elliptic  filters. 

Example  2.6 

Once  again  the  design  example  of  the  previous  section  will  be  used  to 
illustrate  the  new  direct  design  procedure.  The  problem  statement  is: 

Design  an  elliptic  lowpass  filter  that  meets  the  following  specifications: 

•  passband  ripple  of  0.5  dB 

•  passband  ripple-edge  frequency  of  2  kHz 

•  stopband  gain  of  at  most  —20  dB  for  frequencies  greater  than  6  kHz 

•  sampling  frequency  is  20  kHz. 


TABLE  2.9 


ELLIPTIC  LOWPASS  PROTOTYPE  FILTERS 

Passband  ripple  =  0.5  dB  ;  Stopband  gain  =  —20  dB 
(due  to  reference  5) 


N 

0i  p 

02  P 

R 

2 

1.5708 

2.447 

2.7626 

3 

1.5708 

1.9157 

1.9157 

4 

1.5708 

1.6145 

1.0447 

5 

1.5708 

1.5862 

1.0155 

Transfer  Functions 


N 


H(z) 


2 

3 

4 


1 ,575zs  4-2.75  z  +  l  .575 
3.91  zJ4- 1 .1 3z4-l  .22 

1.276z34-2.36  7z24-2.376z4-l  .276 
4.626  z3  +  0 .008  z  2  4- 3 .03  Oz— 0.352 

1.4  179  z4  +2.3  66  z34-3.4362  z7  4-2.336z4-l  .4179 
5.885  z4  —  1 . 0 958  z3  +6. 5858  zJ-l  .0954  z4- 1.338 


1 .794  z8  4-2 .95 12  z*  4-4 .8532  z3 4-4  .85.32  z 


2 


2. 9532*  +  !. 794 


Step  1:  Find  the  critical  digital  design  frequencies. 

2^  =  2*(2  x  10»)  = 

lp  fa  20  x  103  c 

2^,  =  2.(6  xlp  = 

2p  /,  20  x  103 

Step  2:  To  determine  the  filter  order  of  the  prototype  filter,  find  a  and  B2p ,  which 
are  defined  as  follows  from  Table  2.5. 


=  sin  (Bc/2-B'lp/2) 
sin  (0c/2  + 0^/2) 

where  6C  is  the  critical  frequency  of  the  prototype  filter  i.e.,  Bc  =  n  j 2. 
Therefore, 

_  sin(0.57r/2  —  0.2r/2) 
sin  (0.57r/2  +  0.2tt/2) 

sin  (0.7854-  0.31416)  0.45399  (2.49) 

“  sin  (0.7854  -f  0.31416)  “  1.09956 

=  0.5095 

Find  02P,  the  stopband  frequency  for  the  lowpass  prototype  filter. 


= 


j0‘ 

er  *p  —  a 
1  —  ae  2r 

gji.855  _  0.5095 

1  -  (0.5095)e>1855 
1.255e;2'282 


(2.50) 


1.255e~>0-3965 
>B2p  =  2.6785  rad 


=  le>2' 


Step  3:  From  the  design  chart  (Table  2.9),  determine  the  lowest  order  elliptic 
lowpass  filter  prototype.  Table  2.9,  indicates  that  a  second-order  fil¬ 
ter  has  a  stopband  frequency  92p  of  2.447  rad,  which  meets  the  design 
specifications.  Thus,  the  lowpass  prototype  transfer  function  is: 

1.575^2  4-  2.752  4- 1.575 


1 

1 

SteD  4:  Find  the  actual  filter  transfer  function,  using  the  prototype  transfer  func-  ! 

tion  and  the  digital  filter  frequency  transformations  of  Table  2.5.  ! 

HLP(z)  =  HLPp(z)  ; 

z=if^  i 

1.575z2  +  2.75z  +  1.575  (n  * 

3.91.Z2  4-  1.13z  4-  1.22  _  z-o  snsn 

1  —  (  0.  5  09  5  )* 

0.5829z2  +  0.2541z  +  0.5829  j 

3.651z2  -  3.8042z  +  1.6593 

or 

rr  .  .  0.1596z2  +0.0696z  + 0.1596 

Hlp{z)=  z2  -  1.042*  +  0.4545  (2'53)  ! 

The  transfer  function  of  Equation  (2.52)  is  identical  to  the  transfer 

function  obtained  using  the  traditional  design  technique  (see  previous 

1 

section).  ' 

SteD  5:  Verify  the  design  bv  obtaining  a  computer  generated  frequency  response 

plot  (see  Figure  2.16). 

To  further  illustrate  this  direct  design  method,  an  additional  example  J 

, 

follows  involving  the  design  of  a  bandpass  filter.  ] 

Example  2.7 

V, 

A  digital  elliptic  bandpass  filter  is  to  be  designed  to  meet  the  following 

specifications: 

•  0.5  dB  ripple  in  the  frequency  range  600  Hz  <  /  <  900  Hz, 

•  sampling  frequency  of  3  kHz,  and  ! 

is 

ft 

•  maximum  gain  of  -20  dB  for  0  <  /  <  200  Hz.  1 

f 

P 
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Step  1:  Find  the  critical  digital  frequencies  from  the  design  specifications: 

=  =  0.4,  =  1.25  rad 

*  fs  3000 


2tt/u  2tt(900) 


'■  -  T  - = = °-fa  - 188 

^  ^  -  Toir1 = 0133,r  - 0418  rad 

Q'q  =  =  \/(l-25)  (1.88)  =  1.53  rad 

Step  2:  Again,  since  the  design  chart  is  based  on  a  cutoff  frequency  of  9C  =  7t/2, 
the  design  specifications  need  to  be  normalized  to  determine  the  lowpass 
prototype  filter  from  Table  2.9. 

From  the  Table  of  Digital  Frequency  Transformations  (Table  2.5), 
the  following  frequency  transformations  apply: 


Prototype  filter  from  desired  filter: 


s.rm  ^  + tef 

i  - 


(2.53) 


where: 


=  cos(yn/2  +  ^/2) 

“  cos(«i/2-«;/2) 

k  =  tan(0c/2)cot  (9'u/ 2  —  0't/ 2) 

^2P  is  the  stopband  frequency  used  to  determine  the  filter  order,  N. 
9C  is  the  cutoff  frequency  of  the  prototype  filter,  9C  =  tt/2. 


For  this  problem: 


02.  =  ? 


0/2p  =  0.133tt 


=0.6;r 

Of  =0.4tt 


0C  =  0.57T 


COS  (0.37T  +  0.27r) 
COS  (0.37T  —  0.2 7r) 


COs(0.5)7T 

cos(0.l7r) 


a  =  0 


k  =  tan(0.57r/2)cot(0.67r/2  -  0.4tt/2)  =  3.07 
Step  3:  Find  02p  to  determine  the  lowest  order  prototype  filter  that  may  be  used. 


Since  a  =  0 


P-^p  =  _ 


J202  k-l 

_ _  ^  k+ 1 


ej9'p  +  0.509 


1  +  kTiej26'7p  1  +  (0.50 9)e}26'*r 
gjo.836  +  0.509  1.394e”>2-58 

“  ~~  1  +  (0.509)e-'°-836  ~  1.394e>0-282  ~ 

02p  —  2. 86  rad 


=  le->2-86 


(2.55) 


Step  4:  From  Table  2.9,  determine  the  filter  order  of  the  prototype  filter  that 
corresponds  to  the  above  value  of  02p .  It  can  be  seen  from  the  table  that 
a  second-order  filter  fulfills  the  design  specifications. 

Again,  the  lowpass  prototype  transfer  function  is: 


Hlpp{z)  = 


1.575z2  +2.75z  +  1.575 
3.9U2  +  1.13z  +  1.22 


(2.56) 


Step  5:  As  in  previous  design  problems,  the  prototype  filter  transfer  function  is 
converted  to  the  actual  filter  transfer  function  through  the  use  of  the 
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digital  frequency  transformations  of  Table  2.5. 

HBP(z)  =  HLPp(z) 

_  1.575a2  +2.75* +  1.575  (2.57) 

3.91a2  +  1.13a  +  1.22  _  ,n.  m 

-0.509*2  -1 

0.5833a4  -  0.2558a2  +  0.5833 
”  3.6509a4  4-  3.7996a2  +  1.6579 

Step  6:  Obtain  the  frequency  response  plot  for  design  verification  (see  Figure 
2.17). 

In  summary  then,  this  chapter  has  presented  three  commonly  used  recursive 
filter  types:  Butterworth,  Chebyshev  and  elliptic.  The  type  chosen  by  the  designer 
is  dependent  on  the  requirements  of  his/her  particular  design  problem. 

Butterworth  filters  exhibit  a  flat  frequency  response  characteristic,  but  do  not 
have  as  steep  a  transition  region  for  the  same  order  filter  as  do  Chebyshev  and 
elliptic  filters.  With  Chebyshev  and  elliptic  filters  the  steepness  of  the  transition 
region  is  attained  by  accepting  a  degree  of  ripple  in  the  passband. 

All  three  filter  types  may  be  designed  using  either  the  traditional  approach, 
involving  an  analog  prototype  filter  transfer  function,  or  the  direct  design  approach 
that  uses  digital  prototype  filter  transfer  functions. 
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A.  INTRODUCTION 


The  nonrecursive  realization  of  digital  filters  is  desirable  because  it  can  result 
in  two  very  attractive  features:  linear  phase  and  the  absence  of  stability  problems 
because  all  of  the  poles  Eire  at  z  =  0. 

Nonrecursive  filters  are  normally  characterized  by  a  finite  duration  impulse 
response;  consequently,  historical  methods  of  design  involve  the  analytical  deter¬ 
mination  of  the  filter  coefficients  by  expanding  the  desired  frequency  response  in  a 
Fourier  series  and  then  truncating  the  series  to  the  desired  filter  length  (order).  A 
disadvantage  inherent  in  this  design  method  is  the  presence  of  an  overshoot  that 
occurs  at  discontinuities  in  the  desired  frequency  response  (known  as  Gibbs’  phe¬ 
nomenon).  The  use  of  window  functions  was  devised  as  a  remedy  to  this  problem; 
examples  are  the  Kaiser,  Hamming  and  von  Hann  windows  [14]. 

Another  popular  method  of  FIR  filter  design  is  frequency  sampling,  whereby 
the  desired  filter  frequency  response  is  sampled,  and  then  the  inverse  discrete 
Fourier  transform  of  these  sample  values  is  determined  to  find  the  filter  impulse 
response. 

In  this  chapter  both  of  these  methods  will  be  presented,  and  several  design 
examples  given,  but  first  a  review  of  the  theoretical  background  of  nonrecursive 
digital  filter  design  by  Fourier  methods  is  in  order. 
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B.  BACKGROUND 

The  design  of  nonrecursive  digital  filters  hinges  on  the  following  relationships: 
A  causal  nonrecursive  system  can  be  described  by  a  difference  equation  of  the 


form  [1] 


y(n)  =  ^2  b^x(n  ~  k)- 


For  an  input  x(n)  =  ejntf  the  steady  -  state  system  output  is 

yM(n)  =  e>n$H(e>$) 

where  H(eJ0)  is  the  system’s  frequency  response. 


Expanding  the  right  side  of  Equation  (3.1)  for  the  same  input  x(n)  =  ejn  , 
yields  the  following  steady  -  state  output 

y„(n)  =  jr  he*—*’*  =  £ 

Jfc=0  k=  0 

=  bQePn0  +  b^n0e~iB  +  b2ejn9e~>29  +  . . .  +  bLe>nte->Le  (3'3) 

=  e»'n*  [60  +  he~je  +  b2e~j 20  +  . . .  +  bLe~jL0]  . 

Equating  Equations  (3.2)  and  (3.3)  gives 


e>n0H(e]0)  =  e]n0  [b0  +  bxe~i0  +  b2e~j 20  +  . . .  +  bLe~}L e] 


which  implies 


=  Y,  Kt-inS. 


But,  by  definition  the  frequency  response  of  an  LTI  system  is 

n=oo 

H(ej0)  =  J2  Kn)e~in9.  (3.6) 

n=oo 

For  a  causal  filter  ( h(n )  =  0,  n  <  0),  with  a  finite  number  of  delays  (n  =  T), 
however,  the  frequency  response  definition  is 

n  —  L 

*(«*)  =  £  h(n)e-i"'.  (3.7) 


m 

I 

& 

? 

S' 

i 
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Comparing  Equations  (3.5)  and  (3.7)  produces  the  following  relationship  that  is 
the  crux  of  nonrecursive  filter  design 

h{n)  =  bn.  (3.8) 

Equation  (3.8)  establishes  a  direct  relationship  between  the  unit  impulse  response 
of  an  FIR  filter  and  the  coefficients,  £n,  in  the  system  difference  equation. 

The  goal  of  the  design  techniques  presented  in  this  chapter  is  to  determine  the 
filter  coefficients  (or  weights),  b0  ,bx , . . . ,  bi  for  a  desired  frequency  response  charac¬ 
teristic,  Ho{ei0).  These  techniques  include:  Fourier  Coefficient  Design  for  lowpass, 
highpass,  bandpass  and  bandstop  filters;  Windowing;  and  Frequency  Sampling. 

C.  FOURIER  COEFFICIENT  DESIGN 

As  stated  in  the  introduction  we  need  to  determine  the  filter  coefficients,  which, 
in  the  case  of  nonrecursive  filters,  are  also  the  coefficients  of  the  unit  impulse 
response,  h(n).  Thus,  a  relationship  between  the  desired  frequency  response  and 
the  impulse  response  needs  to  be  established. 

Beginning  with  the  desired  frequency  response 

n=  oo 

HD(ej6)  =  £  h(n)e~ine.  (3.9) 

n=  —  oo 

It  can  be  shown,  [1],  that  h(n)  of  Equation  (3.9)  may  be  written, 

2jt+0o 

h(n)  =  (  1/2jt)  I  HD(eje)ejn9dd.  (3.10) 

Je  o 

Equation  (3.9)  is  recognized  as  the  Fourier  series  expansion  of  the  function  H z?(eJ®), 
where  the  h(n )  are  the  Fourier  coefficients  (hence  the  source  of  the  name  for  this 
design  technique). 

Depending  on  the  type  of  filter  that  is  desired,  this  expression  can  be  reduced 
to  the  following  forms: 


V.V.VA'\> 


Lownass  Filters: 


i 

i 

i 

i 

i 


A/,p(n)  =  (/v/tth)  sin(n0c) 

A'  =  desired  magnitude  in  the  passband 
n  =  0,  ±1,  ±2, . . . 


(3.11) 


9C  =  the  desired  cutoff  frequency 

The  number  of  coefficients  is  truncated  to  correspond  to  the  desired  filter  order; 
as  expected,  the  higher  the  order,  the  better  the  approximation  to  the  desired 
frequency  response. 


HD(e 
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Figure  3.1.  Ideal  Lowpass  Filter  Frequency  Response 
Highpass  Filters: 


3 


.'l 

y. 


hHP(n)  =  ^LP(n)(- 1)” 


(3.12) 


Figure  3.2.  Ideal  Highpass  Filter  Frequency  Response 


Bandpass  Filters: 


where, 


= 


*„  = 


^Dp(n)  =  (2  cos(n0o)]  hLP(n) 

eu-ee 


2 

6u  +  ee 


;  0U  =  upper  passband  frequency 
;  6t  =  lower  passband  frequency 
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Figure  3.3.  Ideal  Bandpass  Filter  Frequency  Response 

Bandston  Filters: 

hos{ 0)  =  K  -  hspi 0) 
hBs(n)  =  -hop(n)\  n  =  ±1,  ±2, . . . 


(3.14) 


Figure  3.4.  Ideal  Bandstop  Filter  Frequency  Response 

As  can  be  seen  in  the  previous  impulse  response  expressions  for  the  various 
filter  types,  they  are  all  based  on  a  lowpass  filter  prototype,  h[,p(n).  Thus,  the 
design  steps  for  the  ideal  lowpass  prototype  filter  can  be  summarized  as  follows, 
bearing  in  mind  that  once  the  lowpass  filter  coefficients  are  determined  they  may 
be  transformed  to  the  coefficients  for  a  highpass,  bandpass,  or  bandstop  filter,  if 
desired. 

1.  Lowpass  Prototype  Filter  Design  Procedure 
Step  1:  Translate  the  design  specifications  to  those  of  a  lowpass  prototype.  De¬ 
termine  the  desired  critical  frequency,  9C ,  and  the  passband  magnitude, 
K. 

Step  2:  Find  the  lowpass  filter  coefficients  given  by  Equation  (3.8),  that  is, 

hLp(n)  —  K  =  (A7»n)sin(n0c);  n  =  0,  ±1,  ±2, . . . ,  ±L  (3.15) 

with  L  =  (N  -  l)/2. 

Step  3:  Shift  hip(n)  to  the  right  by  L  terms  to  make  the  filter  causal. 


hLp(n)  =  -  L)Jsin[(n  -  L)tfc];  n  =  0, 1, 2, . . . ,  21  (3.16) 


I 

I. 
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Step  4:  Transform  the  coefficients,  hip(n),  to  lowpass,  highpass,  bandpass  or 
bandstop,  as  desired.  Implement  the  design,  and  compare  the  frequency 
response  obtained  with  the  original  specifications  to  ensure  that  these 
specifications  are  met. 

As  an  example,  suppose  a  highpass  filter  with  a  passband  of  unity  gain  for 
frequencies  greater  than  10  kHz  is  desired.  The  system  sampling  frequency  is  50 


SteD  1:  Cutoff  freauencv.  =  wT  =  =  =  0-4. 

Passband  gain,  K  =  1.0 
Translating  to  lowpass  prototype  specifications: 


=  7r  —  9C  ;  where  9C  is  the  cutoff  frequency  of  the 
lowpass  prototype  and  9'c  is  the  cutoff 
frequency  of  the  highpass  filter 
==>  9C  =  7r  —  0'c  =  7T  —  0.47 r  =  0.6tt 


Step  2:  From  Equation  (3.11),  the  lowpass  prototype  filter  coefficients  are: 

hip{n )  =  sin(0.67rn)/(7rn);  n  =  0,  ±1,  ±2, . . . ,  ±10 

for  a  21-coefficient  filter,  while  the  highpass  filter  coefficients  are  obtained 
from  Equation  (3.12).  The  resulting  frequency  response  is  illustrated  in 
Figure  3.5. 
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n 

hr.p 

Khr 

0 

0.600 

0.600 

±1 

0.303 

-0.303 

±2 

-0.094 

-0.094 

±3 

-0.062 

0.062 

±4 

-0.076 

0.076 

±5 

0.0 

0.0 

±6 

-0.050 

-0.050 

±7 

0.027 

-0.027 

±8 

0.023 

0.023 

±9 

-0.034 

0.034 

±10 

0.0 

0.0 

D.  WINDOWS 

It  has  been  shown  that  the  expression  for  the  desired  frequency  response, 
HD(e>*),  can  be  written  in  terms  of  a  Fourier  series.  Equations  (3.9)  and  (3.10) 
axe  rewritten  below  for  convenience. 


where 


HD(ej0)  =  h(n)e~jn6dd 


2ir+$o 

h(n)  =  (1/2tt)  /  HD(ej9)ejn0d9. 

Je  o 


(3.10) 


In  actual  implementation,  however,  the  infinite  sum  in  Equation  (3.9)  must  neces¬ 
sarily  be  truncated,  since  a  filter  cannot  have  an  infinite  number  of  delays.  Also, 
at  points  of  discontinuity  in  the  ideal  frequency  response  (where  the  magnitude 
abruptly  changes  from  1.0  to  0.0,  or  vice  versa),  the  Fourier  series  approximation 
cannot  match  the  desired  frequency  response  exactly,  even  if  an  infinite  number  of 
terms  were  possible.  There  is  always  an  overshoot  of  about  nine  percent  near 
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discontinuities.  This  overshoot  is  the  well  known  Gibbs’  phenomenon,  and  is  alle¬ 
viated  through  the  use  of  window  functions.  The  filter  weights,  h(n),  axe  modified 
using  one  of  several  available  analytical  expressions  (Hamming,  Von  Hann,  etc.) 
[14]. 


h(n)  =  h(n )  •  w(n) 


(3.17) 


h(n )  =  modified  unit  sample  response 
h(n)  =  original  unit  sample  response 
w(n)  =  window  function 

The  modified  value  of  the  desired  frequency  response  is  therefore 


HD(e>0)=  kn)e-jnB 


n=—  L 


=  ^  h(n )  •  w(n)e  ■7"9 


n=  —  L 


again  where, 


h(n )  =  (1/2 


X)  J  HD(t> 


j6)ejned0 


(3.18) 


which  exhibits  a  gradual  roll-off,  rather  than  the  steep  slope  characteristic  of  the 
ideal  frequency  response.  Thus,  the  tradeoff  involved  when  using  windows  in 
Fourier  design  is  that  a  reduction  in  the  overshoot  caused  by  Gibbs’  phenomenon, 
is  achieved  at  the  expense  of  decreasing  the  sharpness  of  the  frequency  response 
cutoff.  That  is  to  say,  the  passband  and  stopband  ripples  are  suppressed  at  the 
expense  of  a  wider  transition  from  passband  to  stopband;  specifically,  there  is  a 
less  steep  transition.  Figures  3.6  and  3.7  illustrate  these  points;  they  depict  an  un¬ 
windowed  and  windowed  lowpass  filter  design,  respectively.  For  convenience,  some 
of  the  more  popular  window  functions  are  summarized  [14]: 
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Figure  3.6.  Unwindowed  Lowpass  Filter  Frequency  Response 
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Figure  3.7.  Windowed  Lowpass  Filter  Frequency  Response 
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•  Hamming  Window 


w(n)  =  0.54  4-  0.46  cos(mr/L) 
n  =  0,  ±1,  ±2, . . . ,  ±L 


(3.19) 


•  Von  Hann  Window 


w(n)  =  0.50  4-  0.50cos(n7r/L) 
n  =  0,  ±1,  ±2,. . . ,  ±L 


(3.20) 


•  Kaiser  Window 


w(n)  = 


;0(/3[l-(2u/I2)])1/2 


Ud 

n  =  0,  ±1,  ±2, . . . , 


(3.21) 


where  Jo(-)  is  the  modified  zeroth-order  Bessel  function,  and  [3  is  a  constant  that 
specifies  a  frequency  response  tradeoff  between  the  peak  height  of  the  sidelobe 
ripples,  and  the  width  of  the  main  lobe. 


E.  DESIGN  OF  A  DIFFERENTIATOR 

The  Fourier  series  design  procedure  is  easily  applied  to  the  design  of  a  differ¬ 
entiator,  which  is  often  used  in  signal  processing  applications  to  track  the  rate  of 
change  of  a  signal.  The  desired  analog  frequency  response  for  a  differentiator  is 


Hd{jw)  =  jw. 


Its  counterpart  in  the  digital  frequency  domain  is 


(3.22) 


HD(e}°)=j6/T-,0  =  wT. 


Figure  3.8  illustrates  the  ideal  magnitude  and  phase  characteristic. 


(3.23) 


*  k>  v* 


lW 


Figure  3.8.  Ideal  Differentiator  Frequency  Response  Characteristic 
The  Fourier  series  design  approach  requires  the  determination  of  the  Fourier 
series  that  approximates  this  ideal  frequency  response  characteristic. 


Substituting  the  desired  frequency  response  (Eqn.  3.23)  into  Equation  (3.10) 


yields 


h(n)  =  (l/2»)  J  (jO/T)e]n9d6 
=  (j/2zT)  J *  6ejn8de. 


(3.24) 


Evaluating  this  integral,  the  following  expression  is  obtained  for  the  unit  sample 


response  of  a  differentiator. 

h(n)  =  (— l)n/(nT),  n  =  ±1,  ±2, . . . 

h(n)  =  0  ,  n  =  0 


(3.25) 
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Figure  3.9  depicts  the  frequency  response  of  a  20th-order  differentiator  whose  co¬ 
efficients  for  T  —  1  are  as  follows: 
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Also,  included  is  a  20th-order  differentiator  with  a  Hamming  window,  which,  as 
can  be  seen  in  Figure  3.10,  is  a  very  good  approximation. 

In  summary  then,  the  Fourier  series  design  technique  is  very  effective,  and 
can  be  used  in  many  applications  such  as  the  design  of  a  differentiator  illustrated 
here.  The  use  of  window  functions  is  necessary,  however,  to  reduce  the  Gibbs’ 
phenomenon  effect,  which  introduces  a  sacrifice  in  the  designed  filter’s  transition 
region. 

F.  FREQUENCY  SAMPLING 

The  window  design  technique  introduced  in  the  previous  section  has  a  distinct 
disadvantage  in  that  the  computation  of  Fourier  series  coefficients  for  the  desired 
frequency  response  can  be  impossible  when  an  analytical  expression  for  the  desired 
filter’s  frequency  response  is  unknown  or  extremely  difficult  to  determine. 

Frequency  sampling  is  a  technique  originally  proposed  by  Gold  and  Jordan  and 
later  developed  by  Rabiner  et  al.  as  a  solution  to  this  problem  [9].  This  method 
has  two  distinct  advantages: 
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•  It  is  possible  to  design  a  filter  that  approximates  any  desired  frequency  domain 
specifications,  without  ever  having  to  determine  an  analytical  expression  for 
the  filter  frequency  response,  H(e}9). 

•  The  filter  coefficients,  h (n),  can  be  determined  using  the  IDFT  algorithm. 
Before  proceeding  with  an  explicit  summary  of  the  steps  involved,  and  a  de¬ 
tailed  example,  a  discussion  of  the  theoretical  background  will  be  presented. 

In  many  filter  applications  a  sharp  cutoff  amplitude  characteristic  and  linear 
phase  are  desired.  To  ensure  that  these  objectives  are  met  when  designing  a  filter, 
requires  consideration  of  the  number  N  and  location  of  the  equispaced  frequency 
samples.  Discussion  of  these  parameters  relative  to  their  impact  on  the  filter's 
frequency  response  involves  consideration  of  four  separate  cases  [9]. 

Case  1  y  odd,  frequency  samples  at 


9k  =  2xk/N,  fc  =  0, 1, 2, . . . ,  AT  —  1. 


Figure  3.11a.  Frequency  Samples  for  N  Odd 


Case  2  N  even,  frequency  samples  at 


=2xk/N,k  =  0,1,2 . N-l. 
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Figure  3.11b.  Frequency  Samples  for  N  Even 


Case  3  N  odd,  frequency  samples  at 

ek=2T{-~^/-\k  =  0,1,2 . N-l. 


Figure  3.11c.  Frequency  Samples  for  N  Odd 
Case  4  N  even,  frequency  samples  at 


6i.  =  2rr 


(k  -r  1/2) 

N 


.k  =  0,1,2 . .V  -  1. 
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Figure  3. lid.  Frequency  Samples  for  N  Even 
It  should  be  noted  that  Cases  3  and  4  do  not  readily  lend  themselves  to  inverse 
DFT  computation  using  the  FFT  algorithm,  since  the  first  frequency  sample  is  not 
at  0  Hz.  For  this  reason  a  detailed  discussion  of  these  two  cases  will  be  omitted, 
however,  information  regarding  these  cases  can  be  found  in  [9].  A  discussion  of 
frequency  sampling  considerations  for  Cases  1  and  2  follows, 
amnline  Theorem 


Recall  that  for  an  FIR  filter  with  impulse  response,  h(n)  =  {h( 0),  A(l), . . . , 
h(N  —  1)}  ,  the  transfer  function  is 


(3.2C) 


As  discussed  earlier,  the  impulse  response,  h(n),  can  also  be  represented  in  terms 


of  its  Discrete  Fourier  Transform  (DFT),  since  it  is  of  finite  duration  [9] 


h(n)  =  (l/.V)  £  Hke>2*kn'N 


(3.27) 


where  Hk  =  ff(r) 


\z  =  c,1wk/N 


£ 


$ 

S 


Substituting  the  expression  for  h(n)  into  Equation  (3.26)  gives  [9] 


iV-i  r  iV-i 

H(z)  =  ^  (1/JV)  ^2  Hkej2irkn/N\  z~n 

n=0  .  k= 0 

n- i  r/v-i 

=  (1  /N)  Y  Y  Hke>2*kn/N  z~n 


n= 0  L  k  =  0 
N-l  iV— 1 


=  (1  /N)  Y  Y  Hke>2*knlN z~n 

k—Q  n=0 
N- 1  N- 1 

=  (1/N)  YHkY  ej2l'kn/Nz~n- 


k=0  n=0 

Applying  the  finite  geometric  sum  property  to  the  inner  summation 


N —1 

l  - 

ti2rk/N  1 

ej2nkn/N z~n  _ 

z 

l  - 

'  ej  2  rk/N  ' 

n= 0 

z 

N-l 

l  - 

V  ei2"kn/Nz-n  _ 

.  J 

>-[ 

ej2*k/N  ] 

n=0 

*  J 

l  - , 

e}2*kz-fi 

1  _  ej2nk/N  z- 

Substituting  back  into  Equation  (3.28)  yields  [9] 

N-i  „ 


h(z)  =  (i in)  £ 


a ziz!n 

N 


-N\  N~* 


_ fffc _ 

1  _  ej2xk/N 2-l  ' 


Let  2  =  ej9 ,  to  determine  the  interpolated  frequency  response 
...  e~>N9'2  {e^9'2  -  e->n9l2)  ^  Hk 


1  -  (ei2nk>N)  (e~J8) 

rk  ( eiN 9 12  -  e~j9N/2) 

1  _  (ei2*k/N)(e->9) 


a 


e-}N8/2  ffk  (ejW/2  _  e-jN8/ 2)  .  e~j*k/N  .  ejO/ 2 

_  TV  ^  e-jirk/N  .  ej6/ 2  _  ejnk/N  .  e~)8/2 

k=0 

e-jN8/2  .  ej*/2  ffk  (e]N0/2  _  e-;JV9/2j  .  e-jwk/N 

~  jV  ^  ej(8/2-*k/N)  _  eH»/2-vk/N) 

fc=0 

*  fcS  '“(I-#)* 


(3.29) 


Thus,  Equation  (3.29)  expresses  the  interpolated  frequency  response  in  terms  of 
the  sample  values,  Hk ,  of  the  desired  frequency  response  [9]. 

In  Case  1,  where  the  number  of  frequency  samples,  TV,  is  odd,  choosing  the 
frequency  samples,  Hk ,  to  be  real  and  symmetric  yields  a  real  and  symmetric 


impulse  response 


(N-l)/2 


x  Ho  '  2Hk  (2*  \ 
',(n)  =  lv+  I;  —  cos(iT,t/) 


(3.30) 


The  interpolated  frequency  response,  H  (ej8),  is  also  purely  real  which,  as  stated 
in  the  beginning  of  this  section,  is  highly  desirable  for  most  filter  applications 


( N-D/2 

a  («")  =  £  2h ( n ) cos(nd) . 


(3.31) 


n=_  1*^11 


Case  2,  however,  presents  a  problem.  Looking  at  Equation  (3.29),  it  can 
be  seen  that  for  N  even,  the  term  e~jnk!N  inside  the  summation  introduces  an 
imaginary  component  into  the  interpolated  frequency  response,  H  (e}6) . 

For  example,  suppose  the  following  lowpass  filter  frequency  response  is  sampled 
from  — 7r  <  0  <  7r  with  an  even  number  of  frequency  samples,  N  =  4.  From 
Equation  (3.27),  the  following  noncausal  impulse  response  is  obtained 


h(n)  =  (1/4)  £  Hkej2irnk/4 


k=-2 


(3.32) 
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Figure  3.12.  Desired  Lowpass  Filter  Frequency  Response 


therefore, 


h(- 2)  =  (l/4)[0  +  +  1  +  t-’*]  =  -0.25 
k(-l)  =  (l/4)[0  +  e>*2  +  1  +  e~iir/2\  =  0.25 


h(  0)  =  (l/4)[0  +  1  +  1  +  1]  =  0.75 

/i(l)  =  (l/4)[0  +  e~j*/2  +  1  +  e^2]  =  0.25 
and  the  frequency  response  is 

=  £  h( 

n=-2 

=  (~0.25)ej26  +  (0.25)e*  +  (0.75)  +  (0.25)e-^ 

=  (— O.25)[cos(20)  +  j  sin(20)]  +  O.25[cos(0)  -f  j  sin(0)]  +  (0.75)  (3-33) 

+  (O.25)[cos(0)  —  j  sin(0)] 

=  (0.75)  +  (0.5)  cos  8  +  (— O.25)[cos(20)  +  j  sin(20)]. 

The  imaginary  component  introduced  is 

— 0.25j  sin(20). 

According  to  Reference  9,  the  amplitude  of  this  component  should  be, 

a  =  am  £  iw-i)* 

t=- 2 


=  (l/4)[0  -  1  + 1  -  1]  =  -0.25 


which  is  indeed  the  case. 


If  an  odd  number  of  samples,  N  =  5,  were  taken  of  the  same  frequency  re¬ 
sponse,  the  impulse  response  would  be 


h(n)  =  (1/5)  £  Hke»*nk'5 

k=-2 

=  (l/5)(0  +  (lje--'1?1  +  1  +  (lfe-"1*”  +  0] 


(3.34) 


therefore, 


=  (1/5)[1  +  2cos(2ffn/5)] 


h(  0)  =  0.6 


h(±  1)  =  (1/5)[1  +2cos2tt/5]  =  0.324 


h(±  2)  =  (1/5)[1  +  2  cos  4tt/5]  =  0.124. 

The  frequency  response  is  thus 

=  £  h( n)e->* 

n=-2 

=  (0.324)ej2*  -(-  (0.124)ejtf  +  0.6  +  (0.124)e->*  +  (0.324)e~-'2<' 

=  (0.324)[cos(2#)  +  j  sin(2#)]  •+■  (0.124)[cos  9  4-  j  sin#]+  (3-35) 

0.6  +  (0.124)[cos 6  —  j  sin#]  +  (0.324)[cos(2#)  —  j  sin(2#)] 

=  0.6  -|-  (0.648)  cos  9  +  (0.248)  cos(20) 
and  there  is  no  imaginary  component. 

To  remedy  the  imaginary  component  problem  for  the  even  case,  the  following 
substitution  is  made  for  the  frequency  sample  values  Hk  [9] 

Hk  =  Gkejwk'N.  (3.36) 

The  set  of  GVs  is  selected  so  that  G(0)  =  ff(0),  G(N / 2)  =  0, 
and  G{k)  =  -G(N  -*),*  =  1,  2,  . . .,  (JV/2)  -  1. 

The  resultant  impulse  response  is  both  real  and  symmetric,  with  h(n)  =  h 
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(N  —  1  —  n);  n  =  0,  1,  2,  . . .,  (N/2)  —  1.  The  values  of  A(n)  are  given  by, 

/i(n)  =  fr+  Z  ^cos  (^)*+(i7)nfc  (3-37) 


and  the  interpolated  frequency  response  H(ej6)  is, 

( N-D/2 

H(e*e)  =  2h(n)  cos(n0). 

n=_i"=il 


(3.38) 


The  steps  involved,  when  designing  a  filter  using  the  frequency  sampling 
method,  are  as  follows: 

Step  1:  Given  the  continuous  frequency  response  specifications  for  a  desired  filter, 
N  samples  are  taken  at  equispaced  frequencies  over  one  period  of  the 


desired  response 


H(k)  =  H(ej*);d=:2irk/N. 


(3.39) 


If  the  number  of  samples,  N,  is  even,  they  need  to  be  transformed  using 
the  following  relationship,  prior  to  proceeding  with  Step  2 


H(k)  =  G(k)e}27rnk/N 


(3.40) 


such  that 


G(0)  =  H(0) 

G(N/  2)  =  0 

G(k)  =  —G(N  —  k);k  =  l,2,...,  {N/2)  —  1. 

Step  _2:  Determine  the  Fourier  coefficients  of  the  desired  filter  by  finding  the 
ID  FT  of  these  sample  values.  In  other  words,  determine  the  impulse 


response. 


h(n)  =  (1  IN)  H(k)el2""k/N 


(3.41) 
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Step  3:  Again,  it  may  be  necessary  to  use  an  appropriate  window  function, 


h(n)  =  h(n )  •  w(n). 


(3.42) 


Step  4:  The  Fourier  coefficients  thus  determined  correspond  to  the  weighted  filter 
impulse  response,  h(n),  and  the  filter  is  realized  nonrecursively  by  the 
difference  equation 

N- 1 

y(n )  =  ^2  bkX(n  ~  k )  (3.43) 


where 


bk  =  h(n). 


Step  5:  Verification  of  the  filter  design  is  accomplished  by  determining  the  inter¬ 
polated  frequency  response,  ff(eJ®),  resulting  from  the  use  of  the  above 
filter  coefficients.  This  frequency  response  is  compared  to  the  original 
desired  frequency  response,  Hd(^6)i  to  see  if  it  is  a  reasonable  approxi¬ 
mation. 

This  procedure  is  illustrated  by  the  design  of  a  bandpass  filter.  First,  an  odd 
number  of  frequency  samples  will  be  used  and  then  an  even  number. 

Example  3.1  Frequency  Sampling  for  N  Odd 

A  bandpass  filter  is  desired  with  the  ideal  frequency  response  characteristic 
illustrated  below.  The  number  of  frequency  samples  to  be  used  is  N  =  255. 

Sampling  gives 

9  k  —  2nk/N 


where  N  =  255 

k  =  0,1,...,  254 
or  A9  =  2tt/N  =  0.0246. 


(3.44) 
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Figure  3.13.  Desired  Bandpass  Filter  Frequency  Response 


Thus, 

rl.O;  Jfc  =  64,...,  106  and  Jfc  =  149 . 191 

Jf(Jfc)  =  {  0.0;  Jfc  =  0, . . . ,  63  and  Jfc  =  107, . . . ,  148 

l  and  Jfc  =  192, . . . ,  254  . 

Taking  the  IDFT  of  the  above  frequency  sample  values  and  truncating  J'i(n)  to  51 


(  1.0;  it  =  64,... 

=  <  0.0;  Jfc  =  0,...,i 


terms  with  a  rectangular  window,  produces  the  interpolated  frequency  response  of 
Figure  3.14a,  while  Figure  3.14b  illustrates  the  response  obtained  using  a  Hamming 
window. 

Example  3.2  Frequency  Sampling  for  N  Even 

Repeat  the  previous  example,  with  N  =  256  frequency  samples. 

Sampling  gives 


where 


ek  =  2-J t/.v 


N  =  256 


k  =  0,1,...,  255 


A0  =  2-/.V  =  0.0245. 
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Figure  3.14a.  Unwindowed  Bandpass  Filter  Frequency  Response 
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Thus, 


H(k) 


1.0;  Jfe  =  64,...,106  and  k  =  150, . . . ,  192 
0.0;  k  =  0, . . . ,  63  and  k  =  107, . . . ,  149 

and  k  =  193, . . . ,  255  . 


But  recall  for  N  even,  the  sample  values  H(k),  must  be  transformed  to  elimi¬ 
nate  the  unwanted  imaginary  component.  Applying  Equation  (3.40)  gives 


such  that, 


For  example 


H(k)  =  G(k)ej2*k'2!i6 


G(  0)  =  ff(0)  =  0 

G( 256/2)  =  £(128)  =  0 

G(k)  =  —£(256  -  Jfc);  k  =  1, 2, . . . ,  127. 


£(65)  =  77(65)  =  1.0  =  -£(191)  =►  G(191)  =  -1.0 


and  the  transformed  77(65)  and  77(191)  are 

J7(65)  «  i.Oe^65*/256 
77(191)  =  -1.0ei2ir(191)/2M. 

Using  the  transformed  77(7)  values,  the  IDFT  is  determined  to  find  the  impulse 
response,  which  is  again  truncated  to  51  terms.  Figures  3.15a  and  3.15b  show  the 
unwindowed  and  windowed  frequency  responses,  respectively.  It  should  be  noted 
that  when  using  a  Hamming  window  for  an  even  number  of  terms  (in  this  case  52) 
the  equation  for  the  window  function  (Equation  3.19)  should  be  modified  as  follows 


'  (n  _  N-i  \  K 

w(n)  =  0.54  -f  0.46  cos  - - 
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(3.45) 
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Figure  3.15a.  Unwindowed  Bandpass  Filter  Frequency  Response 
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G.  TRANSITION  POINTS 

As  discussed  in  the  previous  section,  frequency  sampling  involves  talcing  N 
equispaced  samples  of  a  desired  filter  frequency  response,  Ho(e^9)-  The  IDFT 
algorithm  is  used  to  approximate  the  unit  sample  response,  h(n),  which,  for  the 
case  of  nonrecursive  filters,  is  equivalent  to  the  Fourier  series  coefficients.  The 
coefficients,  in  turn,  are  used  to  design  a  filter  that  approximates  the  original 
desired  continuous  frequency  response. 

A  filter  designed  in  this  manner  has  a  frequency  response  that  is  equal  to 
the  desired  frequency  response  at  the  frequencies  of  the  sampled  values,  however, 
the  response  can  deviate  significantly  from  the  desired  response  at  frequencies  in 
between  these  sample  values. 

A  method  that  can  be  used  to  smooth  the  frequency  response  involves  the 
use  of  transition  samples  between  the  passband  and  stopband  of  the  desired  fre¬ 
quency  response  [9].  The  values  of  these  transition  samples  jure  selected  based  on 
minimization  of  the  ripple  in  the  passband  and/or  stopband,  or  minimization  of 
the  maximum  sidelobe  of  the  frequency  response.  In  other  words,  a  minimization 
algorithm  is  used  to  find  the  optimum  values  for  the  transition  samples  based  on 
the  selected  minimization  criterion. 

As  before,  an  example  will  be  given  to  illustrate  this  technique,  but  first  a 
summary  of  the  steps  involved: 

Step  1:  Sample  the  desired  continuous  frequency  response  at  N  equispaced  fre¬ 
quencies  where  N  is  the  number  of  Fourier  coefficients  that  will  be  used 
to  approximate  the  desired  filter  response.  The  spacing  between  the 
frequencies  is  AO  —  2n/N. 

Step  2:  Determine  the  impulse  response  h(n)  (Fourier  coefficients)  by  finding  the 
Inverse  Discrete  Fourier  Transform,  IDFT,  of  the  sample  values. 
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Step  3:  Using  the  coefficient  values  determined  in  the  previous  step,  find  the  con¬ 
tinuous  frequency  response  and  compare  with  the  desired  filter  response. 
Step  4:  Use  the  minimization  algorithm  to  adjust  the  transition  coefficient  values, 
in  order  to  obtain  as  close  a  match  as  possible  between  the  desired  filter 
response  and  that  obtained  through  the  design  procedure. 

The  minimization  algorithm  as  used  in  [9],  is  based  on  minimizing  the  maxi¬ 
mum  sidelobe  of  the  frequency  response.  Tables  were  generated  for  lowpass  filters, 
bandpass  filters  and  wide-band  differentiators  of  varying  bandwidths;  1,  2,  or  3 
transition  points;  and  odd  and  even  values  of  the  number  of  frequency  samples,  N. 

\ 

Table  3.1  is  a  reprint  of  subsets  of  these  tables  (due  to  Reference  9)  for  lowpass 
filter  design,  using  1,  2,  or  3  transition  points.  The  number  of  frequency  samples  is 
N  =  15,  the  column  labeled  minimax  refers  to  the  maximum  sidelobe,  and  Type-1 
data  means  that  the  first  frequency  sample  is  taken  at  9  =  0. 

Design  Example 

In  this  example  the  goal  is  to  design  a  lowpass  filter  with  three  frequency 
samples  in  the  passband  symmetric  about  the  origin  (BW  =  3),  and  a  total  of  15 
frequency  samples  ( N  =  15).  The  effects  of  using  0,  1,  2,  and  3  transition  points 
will  be  investigated. 

The  values  of  the  frequency  samples,  transition  point  values,  and  impulse 
responses,  h(n),  axe  summarized  below  for  all  four  cases;  the  transition  point  values 
were  obtained  from  Table  3.1. 
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TABLE  3.1 

LOWPASS  FILTER  DESIGN 
(TYPE-1  DATA,  N  ODD) 


ONE  TRANSITION  COEFFICIENT 

(due  to  reference  [9j) 


BW 

Minimax 

1 

N  =  15 

-42.30932283 

0.43378296 

2 

-41.26299286 

0.41793823 

3 

-41.25333786 

0.41047363 

4 

-41.94907713 

0.40405884 

5 

-44.37124538 

0.39268189 

6 

-56.01416588 

0.35766525 

TWO  TRANSITION  COEFFICIENTS 
BW  Minimax  T\  T2 


1 

2 

3 

4 

5 


N  =  15 

-70.60540585  0.09500122 


-69.26168156 

-69.91973495 

-75.51172256 

-103.46078300 


0.10319824 

0.10083618 

0.08407593 

0.05180206 


0.58995418 

0.59347118 

0.58594327 

0.55715312 

0.49917424 


THREE  TRANSITION  COEFFICIENTS 
BW  Minimax  T\  T2  T3 

N  =  15 

-94.61166191  0.01455078  0.18457882  0.66897613 

-104.99813080  0.01000977  0.17360713  0.65951526 
-114.90719318  0.00873413  0.16397310  0.64711264 
-157.29257584  0.00378799  0.12393963  0.60181154 
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Figure  3.17  depicts  the  results  of  frequency  sampling  design  using  transi¬ 
tion  points.  As  expected,  the  amplitude  of  the  passband  and  stopband  ripple  is 
reduced  with  an  increasing  number  of  transition  points;  however,  the  tradeoff  is 
that  the  sharpness  of  the  cutoff  is  reduced.  This  example  illustrates  the  usefulness 
of  minimization  algorithms  in  smoothing  interpolated  filter  frequency  responses 
obtained  with  the  frequency  sampling  technique.  As  stated  earlier,  Reference  9 
contains  tables  for  lowpass  and  bandpass  filter  design,  as  well  as  wideband  differ¬ 
entiators.  While  quite  extensive,  a  designer  may  want  parameters  that  have  not 
been  tabulated.  In  this  case,  approximate  values  of  the  transition  coefficients  can 
be  derived  from  the  tabulated  values  given,  using  linear  interpolation. 

H.  DESIGN  OF  AN  INTEGRATOR 

In  Section  A  analytical  methods  were  used  to  design  a  differentiator.  The  use  of 
these  methods  to  design  an  integrator,  however,  does  not  yield  an  easily  obtainable 
solution.  As  will  be  shown  shortly,  frequency  sampling  solves  this  dilemma  by 
providing  a  straightforward  design  technique  that  produces  extremely  good  results. 


A  bandpass  integrator  with  the  following  frequency  response  characteristic 
is  desired.  Both  of  the  aforementioned  design  methods  will  be  applied  to  the  design 
problem  to  illustrate  the  advantage  of  the  frequency  sampling  technique. 

a.  Method  1  -  Analytical 

In  the  s-domain  an  integrator  is  described  by  the  transfer  function 
H(s)  =  1/s.  The  frequency  response  is  thus:  H(jw )  =  1/jw  =  (ju>)-1 .  Converting 
the  frequency  response  to  digital  frequencies 


H(jw)  =  1/jw  =»  =  —jT/Q\  T  =  period. 


(3.46) 
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Figure  3.17.  Frequency  Responses  for  Varying  Numbers 
of  Transition  Points 
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Figure  3.18.  Ideal  Integrator  Frequency  Response  Characteristic 
To  determine  the  impulse  response,  h(n), 
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Je>  ~  ' 

-e,  e7 


- 6 , 


h(n)  =  (l/2ff)  f  Hd  (e^)  ej9'c 19  +  1/2*  /  ’  HD  (e")  e* 

Jex 

h(n)  =  (1/2*)  1  (-jT/0)  e^c/0  +  1/2 x  J* 
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ax  er  x“ 
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h(n)  =  -jT/2* 
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A  considerable  number  of  steps  later  yields 

2-2! 


J  o. 


h(n)  =  T/n 


n(d2  -*i)  +  ~^i)  +-  +  “7(92"  -0?) 


n  •  n! 


(3.47) 


where  n  is  the  number  of  coefficients  desired  for  the  implementation  of  the  bandpass 
integrator. 

The  analytical  expression  for  the  impulse  response  of  the  bandpass 
integrator,  Equation  (3.40),  is  cumbersome,  especially  when  the  integrator  design 
is  of  a  high  order.  Therefore,  an  alternative  method  is  desirable. 
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b.  Method  2  -  Frequency  Sampling  ] 

Recall  that  the  expression  for  the  desired  frequency  response  is 


Hd  (e**)  =  -jT/9.  (3.48) 

Letting  T  —  1,  one  hundred  and  one  frequency  samples  will  be  taken  of  the  de¬ 
sired  bandpass  integrator  with  a  bandwidth  of  0.l7r  to  1.67T.  In  other  words,  101 
frequency  samples,  H(k),  of  the  magnitude  and  phase  characteristic  (Figure  3.19) 
will  be  taken  using  the  following  increment  of  the  digital  frequency  0 


AS  =  2tt k/N  =  0.062A:  for  N  =  101 


H(k)  =  H  (eje) 


for  k  =  —50, . . . ,  +50. 

0=0.062* 


(3.49) 


The  101-coefficient  impulse  response  is  obtained  by  determining  the 
IDFT  of  these  frequency  sample  values,  and  a  21-coefficient  rectangular  window  is 
applied  yielding  the  following  bandpass  integrator  coefficients. 


n 

h(n) 

n 

h(ti) 

0 

0.0974 

li 

0.4813 

1 

0.0793 

12 

0.2391 

2 

0.0946 

13 

0.2248 

3 

0.0439 

14 

0.0818 

4 

0.0362 

15 

0.0562 

5 

-0.0562 

16 

-0.0362 

6 

-0.0818 

17 

-0.0439 

7 

-0.2248 

18 

-0.0946 

8 

-0.2391 

19 

-0.0793 

9 

-0.4813 

20 

-0.0974 

10 

0.000 

Finally,  the  interpolated  frequency  response  is  obtained.  As  can  be 
seen  in  Figure  3.20,  the  results  are  good.  (Solid  dots  indicate  the  ideal  frequency 
response.)  Also  included  in  Figure  3.21  are  the  results  obtained  when  a  Hamming 
window  is  used,  while  Figures  3.22  and  3.23  are  the  results  of  using  41  coefficients. 
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This  example  further  demonstrates  the  versatility  and  usefulness  of 
the  frequency  sampling  technique.  It  allows  for  the  design  of  filters  that  may  be 
unattainable  using  traditional  analytic  methods,  and  is  thus  an  indispensable  tool. 
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Figure  3.20.  Unwindowed  21-Coefficient  Bandpass  Integrators 
(Solid  dots  indicate  ideal  frequency  response) 
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Figure  3.21.  Windowed  21-Coefficient  Bandpass  Integrator 
(Solid  dots  indicate  ideal  frequency  response) 
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Figure  3.23.  Windowed  41-Coefficient  Bandpass  Integrator 
(Solid  dots  indicate  ideal  frequency  response) 


The  basis  of  all  computer-aided  design  (CAD)  techniques  is  optimization.  This 
is  to  say,  a  desired  filter  frequency  response  is  approximated  by  a  particular  filter 
whose  coefficients  are  to  be  determined.  The  accuracy  of  the  approximation  is 
evaluated  according  to  some  criterion,  usually  an  error  function,  that  indicates 
how  large  a  disparity  exists  between  the  desired  filter  frequency  response  and  the 
approximated  filter  frequency  response.  Variable  parameters  of  the  approximating 
function  are  then  “adjusted”  to  optimize  the  filter  design  in  terms  of  this  criterion. 

A.  REMEZ  EXCHANGE  ALGORITHM 

An  extensively  used  computer-aided  design  technique  for  linear  phase  FIR 
filters  is  the  Remez  exchange  algorithm.  The  mathematical  basis  for  this  algorithm 
is  the  weighted  Chebyshev  approximation. 

A  summary  of  the  approximating  and  error  functions  for  this  algorithm  follows. 
It  has  been  shown  in  [1]  and  [13],  that  the  frequency  responses  for  the  four  cases 
of  linear  phase  filters  -  i.e.,  even  or  odd  symmetry  with  an  even  or  odd  number  of 
terms  -  can  be  written  in  the  form: 

H  (e39)  =  e-MN-1)/2ej{ir/2)LH  (ej9)  (4.1) 

where  H  (ej9)  is  a  real-valued  function  used  to  approximate  the  desired  filter’s 
magnitude  specifications,  and  the  remaining  terms  approximate  the  desired  phase. 
Table  4.1  gives  the  values  L,  and  the  form  of  H  (e39)  for  all  four  cases. 


TABLE  4.1 

FREQUENCY  RESPONSES  FOR  LINEAR  PHASE  FILTERS 

(due  to  reference  13) 


L 

H  (e>®) 

Case  1  -  N  odd 

(N  —  l)/2 

Symmetrical 

0 

^  a(n)cos(n#) 

n=0 

Impulse  Response 

Case  2  -  N  even 

N/2 

Symmetrical 

0 

Y  b(n)  cos  [0(n  —  1/2)] 

in  — ~  J 

Impulse  Response 

Case  3  -  N  odd 

(N-D/2 

Anti- symmetric 

1 

Y  c(n )  sin(n0) 

n=l 

Impulse  Response 

Case  4  -  N  even 

N/ 2 

Anti-symmetric 

1 

£3  d(n)sin[0(n  —  1/2)] 

fl—  1 

Impulse  Response 

Using  trigonometric  identities,  the  expressions  for  H  (e;*)  in  Table  4.1  can  be 
rewritten  in  the  following  form: 

H  {e}°)  =  Q  ( e j9)  P  (ei$)  (4.2) 

where  Q  (eJ&)  is  a  fixed  function  of  frequency,  9,  and  P  (e^e)  consists  of  a  sum  of 
weighted  cosine  terms  (see  Table  4.2). 
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TABLE  4.2 

FREQUENCY  RESPONSES  FOR  LINEAR  PHASE  FILTERS 


(due  to  reference  13) 


Q  (*') 


P  (ef) 


(N-D/2 

Case  1  1  J2  a(n)cos(n0) 

n=0 

<AT/2)-l_ 

Case  2  cos(#/2)  ^  6(n)cos(n0) 

n=0 

(N-D/2 

Case  3  sin(#)  ^  c(n)cos(n#; 

n=0 

(N/2)-l  . 

Case  4  sin(0/2)  d(n)cos(n#) 

n= 0 


a(n)-\a(n)  =  2/i(^fi-n 


)  for  n  =  1,2,...,  ^ 


6(1)  =  6(0)  +  1/26(1) 

6(n)  =  1/2  6(n  —  1)  +  6(n)  for  n  =  2,  3, . 
6  (iV/2)  =  1/26  (y  —  l) 


£L  -  1 
’  2  1 


c(l)  =  c(0)-l/2c(2) 

c(n)  =  1/2  [c(n  —  1)  —  c(n  +  1)]  for  n  =  2, 3, , 

c(ivfL_l)=:1/2c(^-2) 

c(^i)  =  l/2c(^-l) 


d(n)  = 


d(l)  =  J(0)  -  l/2d(l) 

d(n )  =  1/2  d(n  —  1)  —  d(n)  for  n  =  2, 3, . 

<J(f)  =  l/2d(f  -1) 


«  _  1 
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Since  Q  (ejd)  is  a  function  of  frequency  only,  it  can  be  seen  that  the  approx¬ 
imating  function,  H  (e-7*),  can  only  be  optimized  in  terms  of  P  (e-7*);  that  is,  the 
dependency  of  H  (ejd)  on  the  filter  coefficients  is  contained  in  P  (ejtf).  Thus,  the 
coefficients  of  P  (e-7®)  can  be  varied  to  achieve  an  optimum  filter  design,  and,  the 
true  approximating  function  can  be  generalized  for  all  cases  as, 

L 

P  (e7*)  =  ^  a(n)  cos(nfl)  (4.3) 

n=0 

where  L  =  (N  —  l)/2or  (IV/ 2)  —  1  depending  on  which  case  is  being  considered, 
and  the  a(n )  are  the  weights  from  which  the  filter  coefficients  can  be  determined. 

As  stated  in  the  introduction  to  this  chapter,  a  measure  of  how  well  the  de¬ 
signed  filter  frequency  response  approximates  the  desired  filter  frequency  response 
is  required.  The  weighted  Chebyshev  approximation  uses  an  error  function  defined 
as  follows: 

E(6)  =  W(8)  \hd  ( eje )  -  H  (e-7*)]  (4.4) 

where 

Hd  (ejS)  =  the  desired  frequency  response 
H  (ej8)  =  the  designed  frequency  response 
W(d)  =  weighting  factor 
E{8)  —  error 

In  order  to  see  the  relationship  between  the  filter  coefficients  and  the  error, 
Equation  (4.4)  can  be  rewritten  in  terms  of  the  functions  Q  (e-7*)  and  P  (ej6). 

£(«)  =  W(»)  [Hd  (e**)  -  Q  (e")  P  (e>»)j 

=  mm  («")  -  p  («")  (4'5) 


Letting 


W{8)  =  W{9)Q  (e>*)  and  HD  (e>*)  =  HD  (e>*)  /Q  (e>s) 


MM 


yields 


E(9)  =  W(0)  \hd  (eje)  -  P  (eJ*)]  (4.7) 


Thus,  the  Chebyshev  approximation  that  is  performed  using  the  Remez  exchange 
algorithm  can  be  stated  as  follows: 

Find  the  set  of  filter  coefficients  (determined  from  the  values  of  a(n))  that 
minimizes  the  maximum  absolute  value  of  the  error,  E(9),  over  the  frequency  range 
of  interest. 

E  optimum  =  min 


6  \E(9)\ 


(4.8) 


(coefficients) 

At  this  point,  a  discussion  of  the  weighting  function,  W(9),  is  in  order.  The 
purpose  of  this  weighting  function  is  to  ensure  a  small  tolerance  for  error  in  critical 
frequency  ranges. 

If  W(9)  is  large,  this  means  a  large  deviation  between  the  desired  frequency 
response,  Ho  (eJ*),  and  the  designed  frequency  response,  H  (e-7®),  cannot  be  toler¬ 
ated.  Looking  at  Equation  (4.4)  we  see  that  if  W(9)  is  large,  the  difference  between 
the  desired  and  designed  frequency  responses,  Ho  (ej9)  —  H  (e-7*)  j  has  to  be  small 
to  keep  the  weighted  error  small. 

Conversely,  if  W(9)  is  small,  the  difference  Ho  (ejS)  —  H  (e]6)  can  be  larger 
and  still  meet  the  error  criterion.  Small  values  for  W{9)  would  be  used  in  frequency 
bands  where  close  approximation  to  the  desired  frequency  response  is  not  critical. 

A  copy  of  a  program  employing  the  Remez  exchange  algorithm  [16]  has  been 
included,  and  its  use  is  best  illustrated  with  an  example.  Before  proceeding  with 
the  example,  the  salient  features  of  the  FIR  Linear  Phase  Filter  Design  program 
can  be  summarized  as  follows. 


n'li-.j  »,■  pfc  m  ^  i»*  Ji*  Jr  a* , 
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1.  Main  Program  Functions 

•  Manage  the  Input  - 

NFILT  =  filter  length  in  terms  of  samples,  3  <  NFILT  <  NFMAX 
NFMAX  =£"  128,  but  can  made  greater  or  smaller  by  the  user 
JTYPE  =  filter  type 

1  =  Multi-Passband/Stopband 

2  =  Differentiator 

3  =  Hilbert  Transformer 

EDGE  =  number  of  frequency  bands,  specified  by  9[owtr  and  9upvtr 
(  maximum  allowed  is  10  ) 

FX  =  desired  frequency  response  magnitude  in  each  band,  \Hd  (ej6)  I 
WTX  =  positive  weight  function,  W(9),  in  each  band 
LGRID  =  Grid  Density  (the  grid  density,  default  value  is  16) 

•  Set-up  appropriate  approximation  problem,  based  on  the  desired  frequency 

response,  Hq  and  the  weights,  W(9). 

•  Yield  Output  - 

o  The  coefficients  of  the  best  impulse  response  obtained  from  the  best  cosine 
approximation . 

(’max  "  \ 

min  9  |£(0)|  J 

max 

o  The  extremal  frequencies  where  E{9 )  =  9  |£,(0)| 

2.  Most  Important  Subroutines 
EFF  -  defines  the  desired  frequency  response,  Hd  (e-^) 

WATE  -  defines  the  weight  function,  W(9) 

REMEZ  -  calculates  the  best  Chebyshev  approximation  of  the  desired  frequency 
response 


3.  Important  Note 

It  should  be  noted  that  the  program  output  is  the  impulse  response,  h(n). 
The  user  must  supply  his/her  own  program  to  calculate  and  plot  the  filter  frequency 
response,  which  is  necessary  to  determine  the  suitability  of  the  filter  design. 

4.  Input  Format 

NFILT  JTYPE  #  Bands  LGRID 
EDGE  (Band  Edges) 

FX  (Magnitude  in  Bands) 

WTX  (Weights  in  Bands) 

5.  Bandpass  Filter  Design  Example 

We  wish  to  design  a  bandpass  filter  that  meets  the  following  specifications: 
Passband  with  a  gain  of  one  for  the  frequency  interval  10  kHz  to  15  kHz, 
with  a  system  sampling  frequency  of  50  kHz. 

•  The  digital  frequency  band  is: 

2ir  (ID4) 

9‘  5(10*)  °'4 

2ff(1.5)  (10*) 

)u~  5(10*)  -°-5T 

•  Normalizing  so  that  6  =  ir  corresponds  to  a  normalized  frequency  of  0.5. 

0.2 

ir  0.5 

~  jTT  =>  0u  =  0.3 
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•  Designate  Band  Edges,  Magnitudes,  Weights:  (the  weights  were  arbitrar 
ily  made  equal  for  this  example) 

Band  1 

Edges  -  0.0  to  0.19 
Magnitude  -  0.0 
Weight  -  10.0 

Band  2 

Edges  -  0.2  to  0.3 
Magnitude  -  1.0 
Weight  -  10.0 

Band  3 

Edges  -  0.31  to  0.5 
Magnitude  -  0.0 
Weight  -  10.0 

•  Select  Filter  Length  (  number  of  coefficients  ) 

N  =  21 

6.  Sample  Input  File 
00021  00001  00003  0 
0.0  0.19  0.2  0.3  Band  Edges 

0.31  0.5 

0.0  1.0  0.0 


10.0  10.0  10.0 


Magnitudes 

Weiehts 


The  output  file  below  was  generated  using  the  above  input  file.  Note 
the  output  consists  of  the  filter  impulse  response,  the  input  values  used  and  the 


deviation  between  the  desired  and  approximated  frequency  responses. 


pt'IT’^  T y ° fT T. ^ 17  keb^onsd  frTij 

o  t  o  t  m  \  l  v r i. t ® 3 T 5 m 
ETCHANTE  A  L30e  ITH'I 

n  A  ?l  D  n  A  5  5  felted 

-I'LTT:,',  ’.'’NOT0  =  2 1 


H  *j 
H  (10) 
H  U) 


t*  rn^fJT  c;-^  it*-**-* 

=  li^Q]1!lP  +  AO"=  M  (  -71) 

-  -0.  1037? 33  1^-01  =  ”(  20) 

=  n,  iiaq770^t>fnn  =  ty  )  y  n  \ 

-  -n.  2  31073oi*'-O3  =  vt  18) 

=  -i.  =  q(  1 7 ) 

=  -0. 220  7  1 74 ft?-/}  3  =  Wf  IS) 
*  o.  irso iRsaF+o*  =  M  IS 
=  -0.  1  27  ISRUT^-n  3  =  14) 

=  -o. 22*7002  1  * «.-7  0  =  H(  13 
=  -0.  °S6  IS  no  3?-04  =  Fr  r  1?) 
a  7.  233440227+7n  =  H  (  11) 


LOV1”’  BAND  EDOF 
rjDoco  pi(vd  *-  o  o  v 
O^ST^ED  VAL'TE 
WEIGHTIN'? 

DEVT ATTTV 
DEVIATION  TN  DB 


BAND  1 
0.1007707 
n_  t  q-inr  in 
0.DOOOO0O 

0*  8446704 
-0.7516047 


n\MD  2 

7000000 
7  0  0  7  o  o  n 
,0000000 
1 0  0  n  ^  n  o 
3u.'jft704 
5721715 


BAND  3 
0.  ,3100000 
O .  sn  0  0  00O 
0.0000000 
1  o. OOn^OlO 
0. 3U467H4 
-0.2516047 


rvpcrw^r  ciJc-nnENC'^O — NAY'rv,A  nv  “ruv  oistno  Q’jryv 

0.  0511  36  3’  0 .  10511  31  0.1562480  '  n.ioooono  0.  70700^0 
0.25113^0  O.TD^OOOO  0.310000  0  *.  3443005  0  .  305226  3 
0. 4 4° 2 03  0  O.  5000000 


8.  Results 

Figures  4.1  to  4.3,  illustrate  the  frequency  responses  for  filter  lengths  N  = 
21,41  and  61,  respectively. 

As  expected,  the  frequency  response  improves  with  an  increasing  number 
of  coefficients.  The  cutoff  is  sharper  and  the  magnitude  more  closely  approximates 
0  dB  in  the  passband.  As  with  other  design  methods,  windows  may  also  be  used 
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to  improve  the  frequency  response  characteristic.  Figures  4.4  to  4.6,  show  the  effect 
of  changing  the  band  edges  as  follows,  for  N  =  21  coefficients. 


Figure  4.4  -  Band  1  0.0  to  0.19 

Band  2  0.2  to  0.3 

Band  3  0.31  to  0.5 


Figure  4.5  -  Band  1  0.0  to  0.15 

Band  2  0.2  to  0.3 


Band  3  0.35  to  0.5 


Figure  4.6  -  Band  1  0.0  to  0.1 

Band  2  0.2  to  0.3 

Band  3  0.4  to  0.5 


Figure  4.4  represents  the  original  cutoff  frequencies,  while  Figures  4.5  and 
4.6  show  the  results  of  not  selecting  the  stopband  frequency  cutoff  close  enough  to 
the  passband  frequency  starting  point,  i.e.,  the  passband  is  broadened  beyond  the 
desired  0.2  to  0.3  specifications. 

Thus,  when  selecting  the  stopband  cutoff  frequencies  one  should  choose 
values  as  close  to  the  passband  frequencies  as  possible.  It  is  of  importance  to  note 
that  the  program  does  not  allow  for  the  upper  stopband  frequency  equaling  the 
lower  passband  frequency.  In  other  words  the  following  parameters  would  not  be 


acceptable: 


Band  1  0.0  to  0.2 
Band  2  0.2  to  0.3 
Band  3  0.3  to  0.5 


swot?: 


0.00000 


0.G2832  1.25664  1.88496  2.51328 

FREQUENCY,  THETA 


3.14160 


Figure  4.1.  21-Coefficient  Remez  Bandpass  Filter  Design 


Finally,  Figures  4.7  to  4.9  are  the  result  of  varying  the  weights  of  the 
original  filter  specifications  as  follows  for  N  =  21. 

Figure  4.7  -  Band  1  weight  =  1.0 

Band  2  weight  =  1.0 

Band  3  weight  =  1.0 

Figure  4.8  -  Band  1  weight  =  10.0 

Band  2  weight  =  1.0 

Band  3  weight  =  10.0 


Figure  4.9  -  Band  1  weight  =  10.0 
Band  2  weight  =  10.0 
Band  3  weight  =  10.0 


Again,  as  expected,  the  frequency  bands  with  the  heavier  weight  assigned 
more  closely  approximated  the  desired  frequency  response,  Figure  4.8.  Here  it 
should  be  noted  that  the  weights  are  taken  into  account  relative  to  one  another. 
Figures  4.7  and  4.9  are  identical  because  in  both  cases  the  relationships  between 
the  weights  is  1:1:1. 

In  this  example,  we  have  seen  that  the  Remez  exchange  algorithm  provides 
an  effective  way  to  design  linear  phase  FIR  filters.  However,  a  primary  drawback 
to  this  procedure  is  that  the  CPU  time  required  grows  quite  rapidly  (on  the  order 
of  N)  with  the  filter  order  N. 

For  example,  fifteen  CPU  minutes  are  required  for  the  design  of  a  lowpass 
filter  with  N  =  512  using  a  CDC  6500  [15].  This  is  the  time  required  for  one 
iteration  of  the  procedure.  Typically,  more  than  one  iteration  is  required  to  find 
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technique  to  reduce  the  CPU  time  required  is  desirable.  Reference  15  presents  a 
method  whereby  the  properties  of  a  high-order  filter  can  be  extrapolated  from  a 
lower-order  filter. 


B.  METHOD  FOR  THE  DESIGN  OF  HIGH-ORDER  LINEAR  PHASE 
FIR  FILTERS  BASED  ON  A  LOW-ORDER  PROTOTYPE 

A  detailed  explanation  of  this  method  will  not  be  presented,  however,  a  sum¬ 
mary  of  the  general  idea  behind  this  technique  follows.  The  interested  reader  is 
directed  to  Reference  15  for  details  of  its  application.  It  should  be  noted  that  the 
term  “high-order”  refers  to  filters  with  orders  approaching  N  =  2048. 

The  underlying  basis  for  the  high-order  filter  design  technique  is  the  observa¬ 
tion  that  in  high-order,  multi-passband/stopband  filters,  extremal  frequencies  (i.e., 
ripple  frequencies)  in  the  broad  part  of  these  bands  are  more  evenly  distributed 
than  in  the  transition  regions. 
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Figure  4.10.  Desired  High-Order  Filter 


Figure  4.10  represents  the  desired  frequency  response  of  a  high-order  multi- 
passband/stopband  filter;  and  it  is  noted  that  in  regions  A,  C,  and  E  there  is  a 
uniform  distribution  of  extremal  frequencies,  while  in  regions  B  and  D  the  distri¬ 
bution  is  non-uniform. 
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The  design  of  this  filter  using  a  lower  order  prototype  involves  the  following 


steps: 


•  To  obtain  the  low-order  prototype,  the  uniform  extremal  frequency  regions,  A, 
C  and  E  are  “cut-out”.  The  remaining  non-uniform  regions  B  and  D  are  then 

merged  to  form  the  lower  order  filter  whose  frequency  response  Hd  (ejB  j  is 
shown  in  Figure  4.11. 


HD(ej9) 


Figure  4.11.  Low-Order  Filter  Prototype 

•  The  Remez  exchange  algorithm  is  then  used  to  obtain  the  filter  order,  N,  and 
the  coefficients  required  that  will  meet  the  specifications  of  HD  j .  This 
is  done  by  varying  the  order  and  the  weights,  until  the  order  and  weights  that 
yield  the  desired  passband  ripple  and  stopband  attenuation  are  found. 

•  Once  the  low-order  prototype  filter  is  obtained,  its  extremal  frequencies  are 
plugged  bads  into  the  appropriate  regions  of  the  original  desired  high-order 
filter,  Figure  4.12. 

•  The  extremal  frequencies  are  used  as  initial  values,  and  a  finail  run  of  the 
Remez  Exchange  Algorithm  is  then  made  to  obtain  the  filter  order  and  impulse 
response  of  the  high-order  filter  design.  Again,  the  filter  order  and  weights  are 
varied  until  an  optimum  design  is  obtained. 

In  conclusion  it  should  be  mentioned  that  in  Reference  21  a  program  is  outlined 
that  automatically  designs  high-order  FIR  filters  using  this  procedure.  Approxima¬ 
tions  of  filters  with  orders  up  to  N  =  4096  have  been  achieved  using  this  program. 
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Figure  4.12.  Low-order  filter  prototype  being  placed 
back  into  the  appropriate  regions  of  the 
desired  high-order  filter 

In  conclusion  it  should  be  mentioned  that  in  Reference  21a  program  is  outlined 
that  automatically  designs  high-order  FIR  filters  using  this  procedure.  Approxima¬ 
tions  of  filters  with  orders  up  to  N  =  4096  have  been  achieved  using  this  program. 

This  section  has  illustrated  the  usefulness  and  relative  ease  of  application  of 
Computer- Aided-Design  (CAD)  techniques,  such  as  the  Remez  exchange  algorithm. 

Since  this  algorithm  is  only  applicable  to  linear  phase  FIR  filters,  it  is  worth¬ 
while  at  this  point  to  diverge  to  a  discussion  of  a  popular  algorithm  for  the  design 
of  recursive  HR  filters,  the  Minimum  p- Error  Design  Method. 
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C.  THE  MINIMUM  p-ERROR  DESIGN  METHOD 

Before  proceeding  with  a  discussion  of  this  technique,  it  should  be  mentioned 
that,  in  the  interest  of  brevity,  the  pertinent  equations  and  relationships  used 
are  presented  without  elaboration.  This  is  because  the  mathmatical  basis  for  the 
minimum  p-error  criterion  is  quite  extensive  and  beyond  the  scope  of  this  thesis. 
The  interested  reader  may  find  details  of  its  derivation  in  Reference  19. 

The  minimum  p-error  design  technique  consists  of  a  generalization  of  the  min¬ 
imum  mean-squared  error  design  method,  wherein  the  optimum  filter  coefficients 
are  determined  by  minimizing  the  following  error  function  [18]. 
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Ey(A)  =  [ff(AA)  -  HD(ek)f 


where 


w  (6k)  =  weighting  factor 

H  (A,  9k)  =  approximating  function  (actual  response) 

Hd  (9k  )  =  desired  frequency  response 

9k  =  digital  frequencies  over  the  range  of  interest,  k 
A  =  vector  containing  the  k  independent  parameters 
(i.e.,  the  filter  coefficients) 

Note:  If  the  value  for  p  in  Equation  (4.9)  is  unity,  this  relationship  describes  the 
mean-square  error. 

The  filter  approximation  problem  can  thus  be  stated  as  follows:  Given  an 
amount  of  error  that  can  be  tolerated,  E?p ,  find  the  set  of  k  parameters  (filter 
coefficients)  such  that  the  error  between  the  approximate  frequency  response  and 
the  desired  frequency  response  is  within  the  stated  tolerance. 
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Again,  looking  at  Equation  (4.9)  it  can  be  seen  that  for  large  values  of  p, 
the  approximate  frequency  response,  H(A,0k),  has  to  be  very  close  to  the  desired 
frequency  response,  to  remain  within  the  pre-selected  error  tolerance  value 

of  E2f .  Furthermore,  a  sufficiently  large  value  of  p  results  in  an  optimal  solution 
that  is  very  close  to  the  optimal  Chebyshev  (or  minimax)  solution  [19].  The  method 
of  solution  of  this  approximation  problem,  as  will  be  shown,  is  dependent  on  the 
error  tolerance  value  selected,  E2p ,  the  form  of  the  transfer  function,  H(z),  and 
the  parameter  vector,  A.  First  a  discussion  of  the  transfer  function  is  presented 
because  its  form  has  a  direct  impact  on  several  important  features  in  digital  filter 
design  [18]. 

As  is  well  known,  three  possible  forms  of  the  filter  transfer  function  are  direct, 
cascade  and  parallel  realizations.  With  respect  to  errors  directly  related  to  the 
physical  structure  of  the  filter,  i.e.,  quantization  effects  and  finite  coefficient  size, 
the  use  of  the  direct  form  of  the  transfer  function  in  filter  realization  is  undesirable. 
This  leaves  the  parallel  and  cascade  forms,  with  the  cascade  form  being  selected 
because  the  zeros  of  the  transfer  function  remain  unaltered  in  the  process  of  cas¬ 
cading,  resulting  in  well  defined  stopbands.  Conversion  of  a  filter  transfer  function 
to  the  parallel  form,  on  the  other  hand,  involves  partial  fraction  expansion,  which 
results  in  the  zeros  being  less  well  defined.  Theoretically  this  should  not  be  the 
case,  but,  when  the  parallel  form  is  implemented  digitally,  the  zeros  are  slightly 
different  than  those  of  the  original  direct  form  transfer  function  due  to  the  finite 
precision  of  the  computer. 

Using  the  cascade  form  of  the  filter  transfer  function  has  other  advantages 
including: 

•  stability  tests  of  the  filter  can  be  accomplished  without  lengthy  calculations, 


•  the  frequency  response  is  of  a  simple  functional  form  that  readily  lends  itself 
to  insights  as  to  how  the  poles  and  zeros  (hence,  the  filter  coefficients)  impact 
the  error  function. 

Thus,  the  first  step  in  the  use  of  this  design  procedure  involves  decomposing 
the  proposed  approximating  transfer  function,  H(z),  into  cascaded  second-order 
sections. 


where 


N 


N 


ff(2)=i0n%)=^n 


t=i 


1=1 


z 2  -f  auz  +  a2» 
Z 2  +  b\iZ  +  &2i 


kQ  =  a  positive  normalizing  constant,  and 


(4.10) 


N  =  filter  order  /2. 

The  parameter  vector,  A,  therefore,  consists  of  the  multipliers  used  in  the 
cascade  structure. 


A  =  [an,  «i2,  bu,b12,. . . ,  au,a2i,  bn,  b2i,  ...,k0] 


(4.11) 


The  following  expressions  for  the  frequency  response  and  group  delay  of  the 
filter  incorporate  the  cascade  structure  and  parameter  vector. 

•  Frequency  Response 

N 


TT  C1  +  a2i) cos 6  +  au  +  j(l  -  a2i) sin 8 

“( A’ ")  -  H(e’  )  -  *»  n  +  + 

•  Group  Delay 

r(A,  0)  =  -4 fi(e>»)  =  V  [fie  ( 7 -  2cos0  +  i„+j2sin» 

d8  \(1  +  &2,)cos0  +  bn  -f-;(l  -  b2i)smd 

r  ^  2cos^  +  aii -t-j2sin^ 


(4.12) 


,(1  4-  a2i)cos^  +  aii  +j(  1  -  a2i)sm0j\ 

Equations  (4.12)  and  (4.13)  describe  the  frequency  response  and  group  delay 
of  the  approximating  function,  H(A,0k).  To  determine  if  the  coefficients  of  the 
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A  vector  yield  an  optimum  filter  design  based  on  the  desired  frequency  response 
HD(9k)  and  the  error  that  will  be  tolerated  £2, (A)  involves  finding  the  partied 
derivatives  of  the  frequency  response  and  group  delay  portions  of  the  error  function, 
in  terms  of  the  filter  coefficients,  and  minimizing  them. 

In  other  words,  the  filter  coefficients  comprising  the  A  vector  are  to  be  selected 
so  as  to  minimize  the  following  partial  derivatives: 

3Eg?°|A)  =  E (h)  (a  (A,  St)  -  Hd  (S*))2-'  (4.14) 

9^u~  =  W  ^ 7  (r(A,  «>)  -  /fD(«,))2'-'  (4  15) 

and  similarly  for  dE2pa(A)/dbu  and  dE2PT(A)/dbu,  etc. 

Upon  examining  Equations  (4.12)  and  (4.13),  it  is  that  apparent  these  par¬ 
tial  derivatives  are  not  easily  attainable  because  complex  expressions  are  involved. 
To  remedy  this  the  polar  coordinates  of  the  poles  and  zeros  are  used  for  param¬ 
eters,  rather  than  the  original  rectangular  form.  Thus,  the  parameter  vector  and 
expressions  for  the  amplitude  and  group  delay  are  modified  as  follows, 

•  A  Parameter 


■  tTJ+r  r.-  nj*  vm 


A  [roi ,  $01 1  fpi ,  0p\ , . . . ,  tqi  ,  Ooi ,  Tp,’ ,  Qpi  ,...,fco 


(4.16) 
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rn  r  TT  {1 -2ro,cos(^-0oi)  +  rOi*}1/2  {l-2r0jCOs(^  +  ^0.)  +  r0i*}1/2 

Q!(  A.  (7)  -  rCQ  I  I  «  Icy  I/O 

,t,  {1  -  2r,jCos(S  -  9pi)  +  rpii)  '  {1  -  2rpjcos(9  +  9pi)  +  rpii }  ' 

(4.17) 
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•  Group  Delay 


t(A,8)  =  £  — 


1  -  rpi  cos  (9  -  dpi) 


1  -  rp,  cos  (9  +  8 pi) 


"  L 1  “  2rp« cos  (#  -  0pi)  +  rpi*  1  -  2r>j  cos  (8  +  8pi)  +  rpi3 

_  1  -  roi  cos  (9  -  floj) _ 1  -  roi  cos  (0  +  8oi) 

1  -  2r0j  cos  (8  -  0OI)  +  roi2  1  -  2r0j  cos  (8  +  8oi)  +  roi 2 

resulting  in  the  following  partial  derivatives  of  the  error  function: 

3£l;r(A)  =  EpM**)  s1  «*.«*)  -  hd(«„))2'-' 

ur °%  V' oi 

and  for  dE2pa(A)/d80i  and  dE2pr(A)/d8Qi,  etc.  where 

da  _  r0j  -  cos (8  -  80j) 

droi  1  -  2r0l  cos  (8  -  8oi )  +  roi 3 

roi  ~  COS  (0  +  80j) 

1  -  2rol-  cos  (8  +  8oi)  +  roi3  . 


(4.18) 


(4.19) 


p 

da 

— ■  /V 

r0ism(8-8oi)  ] 

d8oi 

1 -2r0icos(8 -8oi) +roi2  jj 

_ r0jsin(fl  +  fl0j) 

1  -  2r0j  cos  ( 8  +  8oi)  +  roi3 

dr  _  (l  +  rpii)  cos  (8  -  8pi)  -  2rpi  (4.20) 

drpi  (l  -  2rp,  cos {8  -  8pi)  +  rpi j)2 

(l  +  rpj?)  cos  (8  +  flpi)  ~  2rp; 

(l  -  2rpj  cos(0  +  8pi )  +  rp,j)2 
dr  _  rpi(l-rpi2)sm(8-8pi) 
d8oi  (\ -2r pi  cos (8  -  8 pi) +  rpi3)2 

__  rp<(l  -rp<a)sin(0  +  0pi) 

(l  —  2rpi  cos  (0  +  0p<)  +  rpiJ)2 

The  partial  derivatives  da/3rp,,  5a/d0pj,  dr/dr0l,  and  dr/d80i,  are  the  same  as 
the  above  but  with  changed  signs. 

Thus,  the  approximation  problem  amounts  to  the  minimization  of  nonlinear 
functions  (Equations  (4.19)  and  (4.20)),  of  n  variables  (the  poles  and  zeros).  This 
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problem  is  readily  solved  through  the  use  of  the  Fletcher- Powell  algorithm,  the 
details  of  which  can  be  found  in  Reference  20. 

A  program  that  performs  the  synthesis  of  recursive  digital  filters  using  the 
minimum  p-error  criterion  and  the  method  of  Fletcher  and  Powell  for  function 
minimization  is  contained  in  Reference  16.  Included  in  this  reference  are  an  exten¬ 
sive  program  description,  input  requirements,  dimension  restrictions,  and  examples 
to  illustrate  how  the  program  is  used  for  various  types  of  filter  design.  Also  included 
is  a  copy  of  the  actual  program. 

The  reader  is  reminded,  however,  that  to  make  use  of  the  program,  the  transfer 
function  must  be  in  the  cascade  form  mentioned  previously.  The  program  output 
consists  of  an  argument  vector,  X,  and  a  gradient  vector,  Q_,  corresponding  to  the 
minimization  of  the  functions  FI,  F2,  and  F3  which  are  used  by  the  program  for 
the  magnitude  approximation,  the  group  delay  approximation,  and  the  combined 
magnitude  and  group  delay  approximation,  respectively.  Based  on  X,  the  poles, 
zeros  and  coefficients  of  the  cascade  realization  of  the  filter  are  computed.  If  desired 
by  the  user,  the  frequency  response  is  also  given  [16]. 

D.  SUMMARY 

In  summary  then  it  is  appropriate  to  cite  a  recent  paper  by  Little  and  Gowdy 
[17],  that  evaluates  both  the  optimum  FIR  design  method  (that  employs  the  Remez 
exchange  algorithm),  and  the  minimum  p-error  method  used  in  HR  design.  Their 
evaluation  specifically  investigates  convergence  problems  encountered  when  using 
these  iterative  techniques,  and  considers  the  design  of  lowpass,  highpass,  bandpass 
and  bandstop  filters. 


The  following  is  a  brief  summary  of  the  more  important  points  mentioned  in 
this  article,  that  should  be  considered  when  using  these  iterative  techniques. 

1.  Minimum  p-Error  Design  Method 


a.  Advantages: 

•  Extremely  flexible  -  can  be  used  to  design  filters  with  arbitrary  magni¬ 
tude  and/or  phase  characteristics,  as  opposed  to  being  limited  to  lowpass, 
highpass,  bandpass,  or  bandstop  designs. 

b.  Disadvantages: 

•  Nonconvergence  problems  due  to: 

o  a  poor  guess  for  the  initial  parameter  vector  X 
o  finite  wordlength 

•  Uses  large  amounts  of  CPU  (IBM  3081)  time  (up  to  2  minutes  for  higher 
order  filters,  greater  than  N  —  12) 

•  Requires  a  considerable  amount  of  input  to  specify  one  filter 

2.  Optimum  FIR  Filter  Design  Program 


a.  Advantages: 

•  The  program  is  easy  to  use 

•  Consumes  relatively  little  CPU  (IBM  3081)  time,  provided  the  filter  is 
not  of  extremely  high  order. 

b.  Disadvantages: 

•  Restricted  to  the  design  of  the  more  common  filter  types:  multi-band, 
bandpass,  Hilbert  transform,  and  differentiators 

•  Highpass  and  bandpass  filter  designs  were  poor  when  even  filter  lengths 
were  used,  but  good  for  odd  filter  lengths 

•  The  user  must  supply  an  FFT  program  to  obtain  frequency  response  data 
to  verify  the  filter  design.  Convergence  of  the  program  does  not  guarantee 
an  acceptable  design. 


Keeping  these  considerations  in  mind,  both  programs  can  be  used  effectively 
to  design  a  wide  variety  of  FIR  and  HR  filters.  It  is  suggested  that  if  a  designer 
wishes  to  use  either  of  these  programs  this  paper  is  a  valuable  reference  to  assist 
in  identifying  problems  that  may  be  encountered. 


V.  CONCLUSION 
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The  areas  of  recursive  (HR)  and  nonrecursive  (FIR)  filter  design  have  been 
investigated,  while  the  more  predominant  design  methods  have  been  extensively- 
discussed  and  exemplified.  Additionally,  the  prevailing  computer-aided-design  al¬ 
gorithms  (Remez  exchange  and  Fletcher-Powell)  were  also  presented. 

In  conclusion,  points  in  these  three  areas  that  merit  special  emphasis  are 
summarized  below: 


A.  RECURSIVE  FILTER  DESIGN 

•  The  desired  frequency  response  may  be  obtained  using  a  lower  order  filter 
than  if  a  nonrecursive  realization  were  used,  however,  filter  stability  must 
be  considered,  and  a  linear  phase  characteristic  is  not  guaranteed. 

•  The  traditional  design  methods  involving  Butterworth,  Chebyshev  or  el¬ 
liptic  analog  prototypes,  and  the  bilinear  transformation  are  algebraically 
intensive.  The  direct  design  method  presented  eliminates  the  need  for 
determining  an  analog  prototype  and  for  prewarping,  thereby  reducing 
the  number  of  calculations  required. 

B.  NONRECURSIVE  FILTER  DESIGN 

•  Although  a  higher  order  filter  is  required  than  in  recursive  realizations  to 
obtain  the  same  frequency  response;  nonrecursive  filters  are  always  stable 
due  to  the  all  zero  nature  of  their  transfer  functions,  and  can  exhibit  a 
linear  phase  characteristic. 

•  The  filter  coefficients,  h(n),  may  be  determined  analytically  by  expanding 
the  desired  frequency  response  in  a  Fourier  series,  because  the  Fourier 
coefficients  correspond  to  the  filter  coefficients. 


Windows  may  be  used  to  eliminate  the  overshoot  in  the  frequency  re¬ 
sponse  (Gibbs’  phenomenon)  caused  by  truncating  the  number  of  terms 
in  the  Fourier  series  to  a  finite  value.  However,  windowing  decreases  the 
sharpness  of  the  filter’s  cutoff  region. 
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•  Frequency  sampling  can  be  used  when  an  analytical  expression  for  the 
desired  frequency  response  cannot  be  found.  The  IDFT  is  then  used  to 
determine  the  filter’s  impulse  response  from  the  these  sample  values. 

C.  COMPUTER-AIDED  DESIGN 

•  For  FIR  filters,  programs  using  the  Remez  Exchange  algorithm  are  the 
most  popular,  while  for  HR  filters,  the  Fletcher- Powell  algorithm  is  in 
predominate  use. 

•  CAD  is  especially  advantageous  in  the  design  of  extremely  high-order 
filters,  or  filters  with  arbitrary  frequency  response  characteristics. 

As  stated  in  the  introduction,  the  design  methods  presented  here  are  by  no 
means  the  only  ones  available  to  the  filter  design  engineer;  however,  due  to  the 
quantity  of  information  available,  one  is  easily  overwhelmed.  A  need  exists  for  a 
concise  source  of  information  on  the  popular  design  methods  available,  and  how 
they  are  used.  This  has  been  the  intent  of  this  thesis. 
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BAIN  PFCGE AH:  FIE  LINEAR  PHASE  FILIEE  DESIGN  FECGSAB 
AOTHOBS 


JABES  H-  BCCLEL1AN 

DEPAEISENI  OF  ELECTP.ICAI  ENGINEERING  AND  CCBPUTEI  SCIENCE 
BASSACHUSETIS  INSTITUTE  CF  TECHNOLOGY 
CABBS1DGE,  HASS.  0^13i 


T HO£  AS  8.  PAEKS 

DBPAtlBZHT  OF  ELECTRICAL  ENGINEERING 
EIOE  UNIVERSITY 
HOUSTON ,  TEXAS  77001 


LA '■'FENCE  E.  RABINER 

BELL  LABORATORIES 

HOEB AY  HILL,  NEK  JEESEI  07S74 


INPUT: 

NFILT- 


FILTEE  IENGIH 
JTYPE —  TYPE  OF  FIL1EF 

1  =  HULTIPLE  PASSfcAND/STGPHAKE  FIL1EE 

2  =  D IEEE LENT IATCF 

3  =  HILBRB1  1EAKSFC.2E  FILTE2 
NBANDS —  NUHBEE  OF  BANCS 

LGEIC —  GRID  DENSITY,  .ILL  EE  SET  TC  1b  UNLESS 

SPECIFIED  CTHEEKISE  BY  A  POSITIVE  CONSTANT. 


EDGE (2»NBANDS) --  BANDEDGE  ARRAY,  LONER  AND  UPPER  EDGES  FCE  EACH  BAND 
'.'Ufa  A  BAXIfiCE  OF  10  BANDS. 

FX(NBANDS)--  DESIRED  FUNCTION  AEEAY  (CB  DESIEED  SLOPE  IF  A 
DIFFEBENIIATCB)  FOB  EACfc  BAND. 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 


W  T  X  ( N  d A  N  D  S ) —  WEIGHT  FUNCTION  A  EBAY  IN  EACH  BAND. 

0I5FEF.EHT1ATC2,  TEE  KBIGET  FUHCIICo 
PEO POET 10 DAL  TC  F. 


FCE  A 

IS  INVERSELY 


SAHPLE  INPUT  BATA  SETUP: 
32,1,3,0 


32,1.3,0 

0.1,3. 1,0. 2,0. 35 
0.-25, 0.5 
0.0, 1.0, 0.0 
10.  u,  1.0,10.0 

THIS  DATA  SPECIFIES  A  LENGTH  32  BANDPASS  FI LI EL  KITE 
STOPBANDS  0  TC  0.1  AND  C.-2S  TC  0.5,  AND  PASSBAND  F  BOB 
0. ^  TO  0.35  KITH  KEIGHTING  CF  10  IN  THE  SIOPBAuBS  AND  1 
IN  THE  PASSBAND.  THE  GPID  DENSITY  DEFAULTS  TO  1o. 

THIS  IS  IHE  FILTER  IN  FIGURE  10. 


THE  PGLL08ING  INPUT  DATA  SPECIFIES  A  LENGTH  32  FULLEAND 
DIPPSEEKIlAlCt  KITH  SLOPE  1  AND  WEIGHTING  CF  1/1. 

THE  GRID  DENSITY  KILL  EE  SET  TO  20. 

^'2° 

1.0 

1.0 


COBHON  PI2, AD, DEV , I, Z, GRID, DES,KT, ALPHA, IEZT, NEC NS, NGEID 
/COPS/  MIEp,ICUI 


COKBCN  , _ _ _ 

DIBSNSION  IEXT  (oo)  ,AD(66)  ,ALPHA(66) ,X (b6) ,Y (66) 
CIBENSION  H  (66) 

D1BENSION  DtS  (1045)  ,G2ID( 104  5)  , NT  (1345) 
DIBENSIOK  EDGE  (20)  ,FX  (10)  ,KIi  (10)  ,DX?l!lI  (10) 
DOOaLE  PEECISION  fii.PI 
DOUBLE  PRECISION  AD,D£V,X,Y 
DOUBLE  PRECISION  GEB.D 
INTEGER  B£1,fcD2,HE3,dD4 
DATA  BD1,EEI,£CJ,BD*»/1HB,  1HA,  1BN,  1  HD/ 

INF 01=1  IfiACH  (1) 
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C  ICU1=I1BACH  (2) 

P1=4.0*DA1AN  (1.0DQ) 

PI2=2.JDOO*P1 

C 

C  THE  PECGRAfl  IS  SET  UP  FOE  A  BAXIBUB  IEKGIE  CF  128,  BUT 
C  THIS  UPEER  LIB II  CAN  BE  CHANGED  EY  EEDIBENSICNING  THE 
C  ARE AYS  TEXT,  AC,  ALPHA,  X,  Y ,  H  10  BE  NFBAI/2  ♦  2. 

C  THE  ABA  AYS  DES,  GBIC,  AND  W1  BUST  DISENSICNED 
C  lt>  (NFBAX/2  ♦  2)  . 

C 

NFBAX=128 
100  CONTINUE 
JIIrE^O 

PHCGBAE  INPUT  SECTION 

BEAD  («, 110)  NFIX1,JTYPE,NEANES,E3BID 
IF  (NFILT.ES. 0) SICE 
110  FGBB AT  (4 (I3,3X) ) 

IFJNFILI.LE.NFB1I.ANL.NFIIT.GE.3)  GO  TC  115 
CALL  EEECS 
STOP 

115  IF  (NBANDS.LE.O)  NoAN3S  =  1 

GHID  DENSITY  IS  ASSU5ED  TO  BE  16  UNLESS  SPECIFIED 
C THIS VISE 

IF  (LGKID.LE.O)  LGBlD*1o 
JB=2*NBANDS 


BEAD  (4,1.: 0)  (EDGE  (J)  , J=1,JE) 
SEAL i 4 ,120)  FX  (J)  ,J=1,NHANDS) 
BEAD  i  4 , 1201  (WTX  (J)  ,  J=  1 ,  N B  AN  D S) 


HEAD (4 , 120)  (NIX  (J )  , J=  1 , NB aH US) 

120  FGBSAI  (-.fit. 9) 

IF  (JIXPE.G1.0.ANL. JTYPE.LE.3)  GO  TO  125 

CALL  EENCS 

STOP 

125  N£G=1 

IF  (OTYPL.EC.1)  NEG=0 
NO JD=N  F ILT/2 
«UDE=NFILT-2*NCDD 

nfcks=nfilt/2 

IF  (NCDD.  lQ.  I.AttZ.NEG.SC.O)  AiFC»S=A'F  C«S+  1 

SET  UP  THE  CENSE  GEID.  THE  H DBEEB  OF  POINTS  IN  THE  GBID 
IS  (FILIEH  LENGTH  ♦  1)*GEIC  CENSIIY/2 

GBID  (1)=EDGE (1) 

DELF=LGHID*NFCNS 
UELF=C. 5/DELF 
IF  (N EG. EC. 0)  GO  TC  135 
IF  (EDGE  (1) .IT. SELF)  GRID (1) *DELF 
135  CONTINUE 
J=1 
L- 1 

LBAMC=1 

140  FTfsEDGE (141) 

145  TtJP=GBID(J) 

CALCULATE  THE  DESIRED  BAGKITUBE  ESSPCNSE  AND  THE  HEIGHT 
FUNCTION  CN  THE  GRID 

DES  (J) =EFF (I EBP, FX,'«’7X, LB  AND, JTYPE) 

N’T  ( J )  =N A1 E  (1EB£,FX, ATX, LBANC, JTYPE 
J  =  J+  1 

GEID  (J) *TEBF+DELF 
IF  (GRID  (J) .Gl.FUP)  GO  TC  150 
GC  1C  145 
150  GEID  (J-1) =FQP 

DES  (J-1)=EFF  (F UP, FX,KTX,LE AND, JTYPE) 

VI  (J-1) “SATE (PUP, FX, NIX, LfcA HD, JTYPE) 

LBAND=LBAN0*1 
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IFJLEAND.Gl.NBAKES)  GO  1C  160 
SEIC  (J)  =EEG£  (L) 

GC  1C  IhO 
160  NGEID=J-1 

IF  )t>£G.  NE.  NCED)  GC  TC  165 

IF  GEID (NG2ID) -GT.  (0.5-DELE) )  NGuID=NGEID— 1 
165  CCN1INUE 

SET  OP  A  NEW  APP50XIBAIICN  PBCBLEE  WHICH  IS  EQUIVALEN1 
TC  THE  CBI3INAI  EECEIEB 

IF  (NEG)  170.170, 1H0 
170  IF  (HCED.EC. 1)  GC  1C  200 
EC  175  J=1 .NGEIE 
CHANGE=DCCS (?I*GEIt  (J) ) 

DES(J) =DES  (J) /CHANGE 
175  WT (J) =«I JJ) »CnANG£ 

GO  it  200 

180  IF  (NODD.EC.1)  GO  TC  ISO 
DO  165  J= 1, NGEIE 
CH ANGE-DS1N (F1*GEID  (0) ) 

DES (J) =DES  iJI/CHANGE 
185  WT  (J)=WTJJ)  ’•CHANGE 
GC  it  200 

190  DO  1SS  J= 1 , NGEID 

CHAN  GE*DSIi<  (PI2*GSID  (J)  ) 

CES ( J) —DES  (J) /CHANGE 
195  HI (J) =  UT(J) -CHANGE 

IN1IIAL  GUESS  FOE  THE  EXTEZflAl  FEEvUENCIES — ECUALLI 
SPACEE  ALONG  THE  GflE 

200  TEBr-FLOAT (HGBIE-1) /FLOAT (NFCNS) 

DC  2 1 C  J=1, NFCNS 
XI=J-1 

210  IEX1 ( J) =XT*TEHP+ 1 . 0 
IEZ1 (KFCkS+1) =NGKIE 
NE1=NFCKS-1 
N2  =  h£CNS*  1 

CALI  THE  EEBEZ  EXCHANGE  ALGCRITBS  TO  EC  THE  APPECX1BA11CN 
pecbieb 

CALL  5EKEZ 

CALCULATE  THE  IHPULSE  HESFONSE. 

IF  (NEG)  3C0.30C.320 
300  IF  (NCDD.E2.0)  GC  1C  310 
EC  jCS  J=l,Nfll 
NZSJ=N2-J 

305  b  (0) =0.5*AL?BA (NZBJ) 

E (NFCNS)=ALPHA (1) 

GO  1C  3o0 

310  b(1) =0.25*ALPHA (NECNS) 

CO  315  J  =  2 , NS  1 

NZHJ=NZ-J 

NF2J=KFCNS+2-J 

315  H (J) =0.25* (ALPHA  (NZBJ) +AIPHA (NF2J) ) 

HJ NFCNS) =0.5* ALPHA (1) *0.26* ALPHA (2) 

G C  1C  350 

320  IF JNCDD.EC.O)  GC  IC  330 
H  (1) =0. 25* ALP HA (NFCNS) 

H(2i =0.25*ALFHA  « N 3 1 ) 

EC  325  J=3 , Nfl 1 

KZBJ=NZ-J 

HF3J=NFCNS+3-J 

325  H (J) =0. 25* [ALPHA (NZBJ) -ALPHA (NF 30) ) 

H  (NFCNS) *0.5* alpha ( 1 ) -0 . 25 *  ALP H A (3) 
b (NZ) =0. 0 
GC  it  350 

330  H  (1)  =C/.25*ALEaA  (NFCNS) 
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DC  235  J=2,NB1 

NZBJ=NZ-J 

NF2J=NFCNS+2-J 


H (J) =0.25*  (ALPHA (NZ3J) -ALPSA (N  F 2J)  ) 
fa  ( NrCNSj =0.5*ALPHA(1) -0.25* ALPHA (2) 


PBOGBA3  COTPOI  SECIION. 

350  SBT1E(8.360) 

3b0  FCE3AT  (1H1,  70(18*)// 


0  SEITE (8.360) 

0  FCH3AT  (1H1,  70 ( 1H*) // 1 51. 29HFIHITE 
1 13X , 34HLINEAK  PHASE  DIGITAL  FILTEB 
217X,24HEEflEZ  EXCHANGE  ALGCEITHS/) 
IF  (JTTPE.EO.I)  WHITE (6,3oS) 

5  FOBBAI  (22X,  15HBAN0PASS  FIL1EF./) 


IBPUISE  EESPCNSE  (FIE)/ 
DESIGN/ 


365  FOMAI  (221, 15H8AN0PASS  FIL1EF./) 

IF  (JIIPi.EC.2)  WHITE  (8,370) 

370  FOBBAI  (22X,  14HDIFFEF.ENTIATCE/) 

IF  (JIIPE.E'i.  8)  «EI1E  (8,375) 

375  F0E3AI  (20X, 19HHILBEEI  TEANSF CE3EB/) 

WHITE  (8,  j7a)  NF1L1 

378  FOBBAI  (20X. 16HFIIIEB  LENGTH  =  ,13/) 

»'  E 1 1 E  ( 8  ,  J  6  0 1 

380  F0E2AI (15X,28H*****  I3POLSE  RESPONSE  **•**) 

DO  331  J=1,NFCNS 

K=NFILI*1-J 

IF  (N£S.£W. 0)  WHITE  (6,382)  J,B(JJ,K 
IF  (NEG-Sw-1)  ■ HITE  (3, 383j  J,H(J),K 

381  CONTINUE 

382  F0H3AT  (13X,2HH (, 12, 4H)  *  ,E15.8,5H  =  H(»I3,1U>> 
38 J  FCBBA1  (131, 2HH (,I2,4H)  =  ,£15.8,oH  =  -B  (,13,18)) 

IF(«EG.Ei-1.AND.HC£j.SS-1)  SEII t (b , 33  4)  HI 

384  FGBBAT  ( 13X , 2HH  ( ,12 , 8H)  =  C.O) 

CO  450  K=1,N6ANDS,4 

KUP=K+3 

IF (KUP.GT.NBANDS)  KU?=N3ANDS 

WHITE  (8,365)  (oE1.BD2,£D3,£D4,J,J=K,KUF) 

385  F0F8A1  </24X,4  (**A1 ,13.71)  ) 

.HITE  (o,390)  (EDGE (2*J-li  ,J=K,KUP) 

390  FOBflAT  (21, 15HLC.EE  HAND  EDGE.5F14.7) 

■  BITE  (d  ,  395).  (EDGE  (2*J)  ,J=K.K0E)  „ 

395  FOfca.T  (2X,l3BCPEEF  BANE  EEGt,5F14.7) 

IF J JTXPE. NE.  2)  WHITE  (b, 400)  JFX  (J)  ,J=K,KUE) 

400  FCoBAT  (2X,  13HUES1EED  V  AL'J  £  ,2X  ,  5 1  14 . 7 ) 

IF (JTYEE.E0.2l  W£ITE(8,405)  (F  X  ( J)  ,  J  =  K  ,  KOE) 


DO  420  J=K.KUi 


420  CEVIAl (J) =DEV/WII (J) 

HF.ITE  (6,1(25)  (DEVIAT  (J)  ,J=K.KUF) 
425  FOBS AT  UX,SHD£V1ATI0K,6X,5F14.7) 
lFUlYPi.NE.il  GC  TO  450 
EC  430  J=K, KUF 

430  DEVIAT  (J)=x0.0»ALCG1D (DEVIAT (JJ  *FX 
WHITE  (6,435)  (DEVIAT (J) ,J=h.KUP) 
435  FCH8AT  (21, 15HDEVIATI0N  IN  DB.5F14. 
450  CONTINUE 


450  CONTINUE 

DO  45^  J=  1  ,N2 
IJ*IEXI  (J) 

452  GEID (J) =GoID  ( 


GEID (J) =GnID (IX) 

WHITE  (6, **55)  (GRID  (J)  ,J=1,  NZ) 
FOEBA1  (/2X.-.7HEX1EEBAL  FcEOLiLN 
1  (2J.5F12.7l) 

WHITE (6,460[ 

FOEBAT {/IX, 70  (1H*)/1H1) 

GC  TC  100 


NCIES — BAXIflA  OF  THE  EEECF.  CD5VE/ 


460  FOE  BA 
GC  TC 
ENC 


C  FUNCTICN:  EFT 

C  FUNCTION  TO  CALCULATE  THE  DESIF.ED  BAGBIIUDE  RESPONSE 
C  AS  A  FUNCTION  OF  FBECUENCI. 

C  AN  ABBITEAKI  FUNCTION  OF  FF.ECUENCY  CAN  BE 
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APPROXIMATED  IF  THE  USER  EEF1ACIS  TE15  FUNCTION 
WITH  THE  APPROPRIATE  CODE  TO  EVALUATE  THE  IDEAL 
MAGNITUDE.  NOTE  THAT  THE  PARAMETER  FREQ  IS  TEE 
VALUE  Of  NORMALIZED  FREQUENCY  NEEDED  FCK  EVALUATION 


FUNCTION  SFF  (FREQ, FX.WTX, LEAND, JTYPE) 
DIMENSION  F  X  f  5 )  ,  W I X  (5) 

If  (JTYrE. EQ-2)  GO  TO  1 
£FF  =  F  X (LB AND) 

RETURN 

1  EFF-FX (LBAND) *FB£Q 
RETURN 
END 


C  FUNCTION:  NATE 

FUNCTION  TO  CALCULATE  THE  WEIGHT  FUNCTION  AS  A  FUNCIICN 
OF  FREQUENCY.  SIMILAR  TC  IHI  FUNCTION  Eff,  THIS  JUNCTION  CAN 
HE  EEELACEC  HY  A  USER-WRITTEN  ROUTINE  TC  CALCULATE  ANY 
DESIRED  WEIGHTING  FUNCTION. 


FUNCTION  KATE  (FSEQ  ,  FX,  W’T  Y,  IE  AN  E,  JT  YP  £) 


DIMENSION  Fiji)  , »IX  (5) 
IF JJTYFE. EC.2)  GC  TO  1 
UATE-STX (LEAND) 


UATE=S1X (LEAND) 

RETURN 

1  IFJF1  (LBAND) .IT.O.COOI)  GC  TC  2 
UATf-WIX (LEAND) /FREQ 

RETURN 

2  WAIE»B1X  (LBAND) 

RETURN 

END 


SUBROUTINE:  ERROF 

THIS  ROUTINE  WRITES  AN  EHRCB  MESSAGE  IF  AN 
ERBOE  HAS  BEEN  DETECTED  IN  THE  INFUI  DATA. 


SUBROUTINE  ERROR 
COMMON  /CGPS/NIT2B#IOai 
■  SHE  (6,  1) 

1  FORMAT  (44H 
RETURN 


ERROR  IN  INPUT  DATA 


C  SDBBCDTINE:  RZMEZ 

THIS  SUBROUTINE  IMPLEMENTS  THE  REMEZ  EXCHANGE  ALGORITHM 
FGL  THE  WEIGHTED  CHZBYSEEV  AEPRCX IM AT  1C N  OF  A  CONTINUOUS 
FU ACTION  ■ITH  A  SUM  OF  COSINES.  INPUTS  TC  THE  SUBROUTINE 
ARE  A  DENSE  GflID  WHICH  REPLACES  THE  FREQUENCY  AXIS,  THE 
DESIREE  FUNCTION  ON  THIS  GRID,  THE  WEIGHT  FUNCTION  ON  THE 
GRID,  THE  NUMBER  CF  COSINES,  AND  AN  INITIAL  GUESS  OF  THE 
EXTREMA 1  FREQUENCIES.  THE  PROGRAM  MINIMIZES  THE  CBfcYSHEV 
ERROR  EY  DETERMINING  THE  EEST  LOCATION  CF  THE  EXTREMAL 
FREQUENCIES  (P01MIS  OF  MAXIMUM  ERROH)  AND  THEN  CALCULATES 
THE  COEFFICIENTS  CF  THE  BEST  APPROXIMATION. 


SUBROUTINE  EEMEZ 

COMMON  FI2,AD,D£V,X,Z,GSID,DES,W?,ALPHA,IEXT,NFCNS,NGSID 
COMMON  /CO  PS/ NIT I a , TOUT 

DIMENSION  IEXT  <oo).,AD(66)  ,  ALPHA  (66)  .X  (66)  ,Y  (06) 

DIMENSION  DES (104=) ,GEID(1045) ,wl (1045) 

DIMENSION  A  (bb) ,P (b  =  )  ,Q  (o5) 

DOUBLE  PRECISION  FI2, DbUM, ECEK, DTEMP , A, E ,Q 
DOUBLE  PRECISION  DR, DAK 
DOUBLE  PRECISION  AD,DEV,X,Y 
DOUBLE  PRECISION  GEE.D 


THE  P20GBA11  ALLOWS  A  SAXIBUB  HOSEEE  CF  IIEEAIIOWS  OF  25 

ITEflAX=25 
CEVL-~ 1.0 
HZ- AFCNS+ 1 
KZZ  =  NFCNS+2 
KITE»=0 
100  CONTINUE 

IEX1  (K2Z)  =B3EID+1 
tiITEE  =  NI1EH* 1 

If  (NITER. GT.IIBBAX)  GC  TC  HOC 
CC  110  J=1,NZ 
0X1=IEXI (J) 

CTEHf=GHIO (JXT) 

D1EBE=DCCS (01£aP*EI2) 

110  XlJ)=EIEBP 

JET=JNFCNS-1)/15*1 
CO  llO  J=1,NZ 
120  *CIJ)-DIJ,NZ,JSI) 

CHuB-C. J 
EE£N=Q.  0 
K-  1 

EC  130  J=1,NZ 
1-1121 (J) 

ET£SE=AD (J) *D£S !L) 

DNOB =0  Ai  OB  +E1£BP 
CT£BE=FLCA1 (6) *AD (J)/ST (!) 

CEEN*CE£ft*01£B? 

130  K*-K 
EEV=DHUS/DD£N 
WHITE  (8,131)  DEV 

131  FGEBAT  <  IX,  12HDEVIAIIUE  =  , F12.9) 

IF  "(LEV. GT. 0.0)  N0*-1 
EEV=-FLCA1  (NO) *E£V 
K  =  KU 

CO  14C  J=1,NZ 
L=IE21  (J) 

CI£«E=FLOAT  (K) *EEV/»T  (I) 

X ( J) =E£S  II) *CjlBr 
140  F.i-k 

IF.(DEV.GT.DEVl)  GC  IC  150 
CALL  OUCE 
GC  IC  400 
150  EEVl=OtV 
0CHNGE=0 
F  1=I£XI J 1 ) 

KNZ=IEX1 (HZ) 

KLCW=0 
N01=-N0 
J-  1 


SEAECH  FOE  THE  EXTSEflAL  FEECUFNCIES  CF  IHE  BEST 
A?  EEC  XI  BAUCH 

200  IF  (J.EO.NZZ)  YNZ=CCBP 

t x  f  h’?!  i  t r.  'n'\n 


-NZZ) 

IFlJ.Gt.NZZj  GU  TC  300 
K0E=IEXT (J*1) 


L=I£X1 (J) *1 
NUI=-NUT 


IF  (J.EC.2)  Y1=CCBE 
COEP=DEV 


COBP=DEV 

IF  (L.GE.KOP)  GC  TC  220 
E3E  =  GEE  (L.NZ) 

E8£- (EHc-DEJ  111) *K T  (L) 
CTEEr  =  FLCAl  (NUi  *EfiE-CCaP 
IF  (CTEBP.LE.0.0)  GC  1C  220 
CGBP*FLGAT  (NUT) *EBB 
210  1=1*1 

IF  (L.GE.KUP)  GO  TC  215 
EBE  =  GEE  (L,N2) 


i.k*.  u  ‘  Jk J 
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ESE  =  (EHF-DES(L)  )  *SI  (I) 
£IESF=5L0A1  (SOI)  *Eor-CC3P 
IF (EIEBP.I.E.0.0)  GC  TC  215 
C0EP=F10AT  (SOI)  *£Bfi 
GC  1C  210 
I  Ell <J)=i-1 
J=J  +  1 
K1CS=I-1 
JChKGE=JCHKGE+1 
GC  IC  200 
1*1-1 
1  =  1-1 

IF (I.IE.KLCK)  GC  1C  250 
££F*G££ (1,  NZ) 

£Sr  = (ERF-DES (1) ) *MI (L) 
CT£HF=FL0A1  (KOI)  »Efiii-CCHP 
IF  (EIESF.G1.0.0)  GC  IC  230 
IF  (JCBSGE.1E.0)  GC  TC  225 
GC  1C  260 

CCf1P=F10Al  (SOI)  *EEE 

1=1-1 

IF  (L.1E.KLCW)  GC  TC  2«0 
EBP =G££ (i,K2) 

EBfi= (Ebc-DES (I)) *KT  (1) 

CTEfl  P=FLCA1 (Nil.) *E£P-CCHP 
IF (EIEHP.LE.G.O)  GC  TC  2U0 
CUKP=F1CAT  (SOI) *EEE 
GC  1C  235 
KlCk=l£IT (J) 

IEX1 (J) =L*1 
J=J*1 

JCHSGE=JCHKGE+ 1 
GC  1C  200 
1=IEZ1  (J) *1 

IF  (JCEuGE.GT.O)  GC  TC  ^.15 

1  =  1+1 

IT  (L.GE.KOP)  GC  TC  2oO 

EI.E=OEE  (L.KZ) 

El(f=  (EEi-DES  (1)1  *K1  (1) 
£TinP=FLCAT  (1.01)  *£ER— CCHP 
if  iciEar.ii.c.o)  gc  lo  255 
CCfiP=I 1C AT  (NUI)»EEE 
GC  1C  210 
K1C'»=I£X1  (J) 

j=0  +  1 

GC  1C  200 

If  (J.Gi.HZZ)  GC  TC  520 
IF  (K1.GI.IZXT  MI)  K1=I£XT(1) 

IF  1knZ.LI.IEX1 (NZ) )  KKZ=I£XT  (El 
SOX  1  =  SOT 
MJT=— KU 
1=0 

KUF*  K 1 

COHP=XMZ* (1.00001) 

IOCK= 1 

1=1*1 

IF(l.GE.KOE)  GC  TC  315 

EEF=GEEJL,SZ) 

lui- (FER-DES (1) ) *WT  (1) 

CT£BF=  FICAT ( NO  I ) * E Bt-CCflP 
IF (DTZHP.LE. 0.0  GC  IC  310 
CCnF=F10Al  (NOT) *£s£ 

J=S2Z 
GC  TC  210 
LUCK  =6 
GC  1C  325 

IF  (LOCK. GI. 9)  GC  TC  350 

IF  (CCKP.Gl.il)  I1  =  COflP 

K  1  =  IEil (NZ2) 

l=NGEIO*l 

KIC  *  =KS  Z 

SOT  =  -S  UT 1 
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coap«n»  (i.oojoi) 

330  1=1-1  .  „ 

IF  (1  .IE. KICK)  3C  1C  J40 
E5MGEE  (L.NZ) 

EEF=  (ESF.-EES  (L)  )  *H1  ,'L) 

DTEH t=FLCAT (NUIl •EBB-CCfiP 
IF  (21EHP.it. 3. 0)  GC  IC  330 
U  =  NZ2 

C0Hi=FI0Al  (NOT) *E3B 
10CK=IUCK+  10 
SC  IC  235 

340  IF  (LUCK . EC. 6 1  SC  TO  370 
CO  345  J=l, NFCNS 
NZ2HJ=NZZ-J 
NZEJ*NZ-J 

345  IEXT  (NZZHJ) =  IEX1  (NZHJ) 

1EX1  ll) =K1 
GO  IC  1 00 

350  K  N  =  IE JI  (NZZ) 

DO  360  J=1, NFCNS 

360  IEX1 iJ) *IEil (J  +  1) 

IEXT  (N2) =Kk 
3C  IC  100 

3 70  If  (JCB»G£.GT.0)  GC  1C  100 

CALCULATION  OF  THE  COEFFICIENTS  Of  THE  EE3I  APPECXIBAIICN 
USING  THE  INVE3SE  CISCEZ1E  FCOSIEE  TBANSFCiB 

400  CONTINUE 

Nfl  1  =  NF  CNS- 1 
FSH=  1. 0E-06 
G1EB  F  =  GSj.D  JI) 

1  (NZZ) =-2.0 
CN=2  * NFCNS- 1 
CELF= 1 . 0/CN 
1=1 
KKK=0 

If  (OF ID  (1 )  .LI.  0.  0  1 .  AND.  GF.IC  ( NGEI3)  .GI.0.49)  KKK  =  1 
IF  (NFCNS. LE. 3)  K K K  = 1 


SB*-  (CIEHP  +  CNOaj/ (CIEBP-DN0H) 
405  CONTINUE 

DC  430  J=1, NFCNS 
FT=J-1 
fl=FT=DELF 
XI=CCCS iFI2*Fl) 

IF  (KKK.SS. 1)  GO  IC  410 
Xl  =  (IT-BEJ/AA 
XT1  =  SCET  (  i  . 0-XT* XT) 
FT=A1AN2UI1,XI)/FI2 
410  XE  =  X  (L) 

IF  (ZI.GT.XE)  GC  IC  420 
If  (  (XE-XI) .11. FSB)  GO  10  4  15 
1  =  1+1 
GC  IC  410 

415 

420  IF  IJII-XE) .II.FSH)  GO  1C  415 
GFID  < 1 ) =F1 
A(J)=GEE(1.N2) 

425  CONIINOE 

If(L.GT.I)  1  =  1-1 
430  CONIINOE 

GEID ( 1) =GTEBP 

EDEN  =  PI 2/C  N 

CO  510  J=1, NFCNS 

E1EHF=0.0 

cnub=j- l 

£t.Un  =  DNUa*LOEN 


I 

i 


\ 
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If (Nnl.LT.1)  3C  TC  505 
EC  500  K=l,isl 
EAK  =  A <K  +  1) 

EK  =  K 

DT£KE=DT£BE+DAK*ECCS  (DNUH*EX) 

CIE8E  =  2.0*£lEaP-fA  (1) 

ALPHA iJ) =EIE8P 

EC  550  J=2,HFCNS 

ALPHA (J) =  2 . 0* ALPHA  !  J) /CS 

ALEHA ll) *A1PBA (1)/CN 

If IKKK.EC- 1)  30  IC  545 

P  J1)=2.0»ALeHA(NFCNS)  *Ba* ALPHA  (Hal) 

? 12)  =2. 3* AA* ALPHA (NFCNS) 

0(U=aI-PHA  (NFCNS-2) -ALPHA (SFCNS) 

EC  540  0=2, uni 

If  (J.1I.BB))  GO  1C  515 

AA=0.5*AA 

E3=0.5»BB 

CC Nil HUE 

t  (j+ii=a.o 

EC  =20  K=1,J 
A  (X)  =F  <K) 

P (K)=2.0»EE*A(K) 

EC  525  K= 1 ,  JH  1 

E(K)=P  jS)+2<K)*AA*A(K+1) 

OP  1=J+  1 

EC  530  K=3,JP1 
P (K) =E (K)  +AA*A  (K-1) 

If  <J.EO-  Nfil)  30  TC  540 
EC  535  8=1, J 

Sjli=NFiNa-1-J 

ttlk 

DC  543  J=1,N?CNS 
ALEHA  I  J) =E <J) 

CONTINUE 

IF  (HICSS.G1.3)  EfIUBJ. 

ALPH*  (NFCNS*1)=C.0 
ALEHA  NFCNS+2) =C-0 
2EIUBN 
ENE 


FUNCTION:  D 

FUNCTION  TC  CALCULATE  THE  LAGEANGE  INTESFCLATICN 
COEFFICIENTS  FOB  USE  I»  THE  FUNCTION  GEE. 


DOUBLE  PRECISION  FUNCTION  E(K,N,fl) 

CCnnCN  PI2, AE,DEV,X,  X,GKIE,DE3,.'T,ALFUA,I£XT, 
EI8ENSI0N  iiXIjbo  ,A2(bo)  , ALPHA  INN)  .X  (bb)  „  X  (b 
CI8E  N5ICN  EES (1043) , 5EI2 < 1 C45)  ,21  (1045) 

DOUBLE  PRECISION  AE,D£H,X,X 
DOUELE  PBECISiON  C 
DOUBLE  PRECISION  EI2 
C=1.0 
C=X JK} 

DC  j  1=1,8 
CO  2  J=L,N,8 

Te  j  i  _im  i  r  i 


[  ,  NFCNS , HGEID 
(b6) 


IF  jJ-K)  1,2,1 

1  D=  =  . C»D« (w-X <JJ  ) 

2  CONTINUE 


3  CONTINUE 
D=1. 0/C 
EETUEN 
END 


FQNCIICN:  GEE 


FUNCTION  TC  EV AIUAT I  IBE  FRECOENCT  RESPONSE  USING  THE 
LAGRANGE  INTERPOLATION  FCBMUIA  IK  THE  EAEXCENTEIC  ECBfl 


DOUBLE  PRECISION  FUNCTION  GEE(K.N) 

COMMON  PI2.AE,D£V,X, Y,GRIC,EES,  ALPHA, IE2T , NFCNS, NGF 

DIMENSION  ISXI  (fco)  ,AL(Ho)  , ALrHn (bo)  .XJofc) ,X  (oi) 
ElflEKSICK  DES i iOUS) -GRID ( *G45) ,  iil  (1045) 

ECU3LE  PEECISIOK  P.C.E.Xf 
DOUBLE  PRECISION  PI2 
DOUBLE  PEECISIOK  AD,DEV,X,X 
P  =  C.O 
XI*GEID  (K) 

XF=DCCS  jpl2*XF) 

E=C.  0 

DC  1  J  =  1  ,  H 
C=Xr-X  (J) 

C  =  AD  (J)/C 
E=E*C 

1  r=P*C*Y{J) 

GEi=P/D 

RETURN 

END 

SCBECUTIKE:  CUCH 

•  BITES  Afi  EER OF.  MESSAGE  KHEK  THE  ALGORITHM  FAILS  IC 
CONVERGE.  THERE  SEEM  TO  BE  T.C  CCHEIIICA'S  UNDER  AfilCH  ' 
TEE  ALGO EXT  HE  FAILS  10  CONVERGE:  (1)  THE  INITIAL 
GUESS  FOR  THE  EXTREMAL  FRECOENCIES  IS  SC  PCCR  THAT 
TEE  EXCHANGE  ITERATION  CANNOT  GET  STARTED.  CP. 

(2)  NEAB  THE  IESMINATICN  OF  A  CORRECT  DESIGN, 

THE  DEVIATION  EECBEASES  DUE  TC  BOUNCING  ERRORS 
AND  THE  PROGRAM  SICES.  IN  THIS  lATIER  CASE  THE 
FILIEB  E2SIGN  IS  PROBABLY  ACCEPIAEIE,  EUT  SHOULD 
BE  CHECKED  BE  COMPUTING  A  f EEC  DENCY  EESRCNSE. 

SUBROUTINE  CUCH 


COMMON  /OCRS/ NITER, ICD1 
•  RITE  (6,  1)  NITER 

1  lOLSAI  H^B  »***•»*»**•*  FAILURE 


l  lOLSAI  H4E  ••••»*••**•»  FAILURE  TC  CONVERGE  *•»»»*»***/ 

1 4  IHUrRQBAELE  CAUSE  IS  MACHINE  ROUNDING  ERROR/ 

223H0NUMBER  CE  I1EEATICNS  =,14/ 

339H0IF  THE  NUMBER  OF  ITERATIONS  EXCEEDS  3,/ 

4 o2B0XHE  DESIGN  SAX  BE  CORRECT,  BUT  SHCU1E  BE  VERIFIED  KITH  AN  FF' 
RETURN 
END 
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